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A 

area,  m 

a 

inner  radius  of  vapor  space,  m  (Section  I) 
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coefficient  in  the  discretization  equation  (6.9)  or  coefficient  in 

eqn.  (5.19) 
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term  in  eqn.  (5.12) 
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term  in  eqn.  (5.15) 

b 

inner  radius  of  heat  pipe,  m  (Section  I) 
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term  in  eqn.  (5.19) 
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source  term  in  eqn.  (6.9) 
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Bi 

V/H 

B2 

L/H 

c 

C°/c( 

c° 

coefficient  in  eqns.  (3.5  and  5.9),  J / (kg  -  K) 

Cs* 

cJci 

c 

specific  heat,  J/(kg-K) 

c 

outer  radius  of  heat  pipe,  m  (Section  I) 

cm 

specific  heat  of  mushy  phase,  l/2(cg  +  c^),  J/(kg  -  K) 

CP 

specific  heat  at  constant  pressure,  J/(kg-K) 
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c 

P 

equivalent  specific  heat,  J/(kg  -  K) 

Cv 

specific  heat  at  constant  volume,  J/kg-K 

» 

D 

inside  diameter  of  the  circular  pipe,  or  diagonal  distance,  m  (Section 

III) 

f 

D 

diameter  of  the  vapor  space,  m  (Section  IV) 

D 

diagonal  distance,  m,  also  used  as  dimensionless  conductances  in  eqn. 

(5.19)  (Section  V) 

D 
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diameter  of  the  vapor  space,  m 
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E 

Et 

erf 

F 
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6 

St 

H 


H 

H 

H 

h 

h 

h 

K 

K 

K 


si 


lfg 


2  2 

total  energy  of  the  vapor  per  unit  volume,  p  [cyTv  +  1/2  (U  +  Vy)j 
(Section  IV) 

enthalpy,  J/kg  (Section  V) 

o 

total  energy  of  the  vapor  per  unit  volume,  p(CvT  +  1/2  U  ) 
error  function 

dimensionless  flow  rate  through  a  control- volume  face 

2 

friction  coefficient  at  the  wall,  r/pU 

friction  coefficient  at  the  exit  of  the  evaporator 

friction  coefficient  for  the  fully  developed  turbulent  flow 

2 

gravitational  acceleration,  m/s 
integral  in  eqn.  (3.18) 
integral  in  eqn.  (3.19) 

height  of  thermal  cavities  or  height  of  the  liquid  level,  m  (Section 

II) 

latent  heat,  J/kg  (Section  III) 

height  of  the  vertical  wall,  m  (Section  V) 

3 

latent  heat  of  melting  per  unit  volume,  J/m 
enthalpy,  J/kg 

heat  transfer  coefficient,  V/m  -  K 
latent  heat  of  evaporation,  J/kg 
permeability, 

dimensionless  thermal  conductivity,  k/k^  (Sections  III  and  V) 

Knudsen  number 

Vk« 

thermal  conductivity,  W/(m-K) 
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Nomenclature  (continued) 


L  half  width  of  the  heating  element,  m  (Section  II) 

L  length  of  the  pipe,  m  (Section  III) 

L  length  of  the  heat  pipe,  m  (Sections  IV,  VI,  and  VII) 

L  latent  heat,  J/kg  (Section  V) 

L  length  of  the  adiabatic  section,  m 

a 

Lc  length  of  the  condenser,  ra 
Lg  length  of  the  evaporator,  m 
l  reference  length  in  Figs.  3.2  and  5.5,  m 
M  Mach  number,  U/V7R.T 
M  molecular  weight,  kg/Kmol  (Section  IV) 

Ma  Mach  number,  w/V?p /pv 
m  total  mass  of  the  PCM,  kg 
m  mass  flux,  kg/m  -s 

2 

m^  mass  flux  at  the  liquid- vapor  interface,  kg/(m  -2) 

2 

mQ  rate  of  evaporation  or  condensation  per  unit  area,  kg/m  -s 
Nu  Nusselt  number,  qH/k(Tji  -  Tc) 
n  unit  outward  normal  direction  (Section  IV) 
n  normal  direction  of  the  liquid- vapor  interface  (Section  VI) 

n  time  step  (Section  VII) 

P  (p  +  pQgy)  $ lp<?  (Section  II) 

P  dimensionless  pressure,  (P'P0)//°fU0  (Section  III) 

P  (p  +  /9gy)  H 2 /p  a2],  (Section  V) 

2 

P£  reference  pressure  for  the  Clausius- Clapeyron  equation,  N/m 

2 

Pcr  reference  pressure  for  the  Clausius- Clapeyron  relationship,  N/m 

2 

Pg  pressure  at  the  evaporator  end  cap,  N/m 
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Nomenclature  (continued) 


2 

Pgk  sink  pressure,  N/m 

9 

Pgo  source  pressure,  N/m 

Pr  Prandtl  number,  v/a 
2 

p  pressure,  N/m 

•  * 

p  modified  pressure 

p'  pressure  correction 

*  pQ  inlet  pressure,  N/ni  (Section  III) 

pQ  centerline  pressure  at  the  end  cap  of  the  evaporator,  or  reference 
pressure,  N/m  (Section  VI) 

Q  heat  energy,  J  (Section  III) 

Q  total  heat  input  rate  at  the  outer  pipe  wall  at  the  evaporator,  V 
(Section  VI) 

Q  heat  input,  V  (Section  VII) 

Qc  maj[  maximum  capillary  heat  transport  rate,  ¥ 

(Qr)  heat  transport  factor,  ¥-m 

total  latent  energy  stored,  J 
Qm  energy  storage  density,  J/kg 

* 

Qt  total  energy  stored,  J 
q  heat  flux,  V/m^ 

2 

1  q'  new  heat  flux,  V/m 

2 

q"  new  guessed  heat  flux,  V/m 

2 

q^  latent  heat  flux,  V/m 

3 

qg  heat  source,  V/m 
R  gas  constant,  J/kg-K 

Rq  outer  pipe  wall  radius,  m 
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Nomenclature  (continued) 


Ru  universal  gas  constant,  KJ/Kmol-K 

Ry  vapor  space  radius,  m 

Ry  wick- wall  interface  radius,  m 

Ra  Raleigh  number,  g/TH^T^  -  Tc)/i/a  (Section  II) 

Ra  Raleigh  number,  g/2H^(T£  -  T°)/i^a^  (Section  V) 

Raj  Raleigh  number, 

Re  axial  Reynolds  number,  /?UDv//i 

Re^  fluid  Reynolds  number,  UQ  D/i/^ 

ReQ  radial  Reynolds  number  at  the  wall,  /?V0Dy//i 
r  radial  coordinate,  m 

rc  effective  capillary  radius,  m 

rp  radius  of  curvature  of  the  meniscus,  m 
s  S°/C<  (T?„-T°)  (Section  III) 

S  S°/c^  (T°  -  T°)  (Section  V) 

S  interface  position  in  vector,  m 
S°  term  in  eqns.  (3.5  and  5.9),  J/kg 
St  Stefan  number,  c^(T°n-T°)/H  (Section  III) 

St  Stefan  number,  c^(T^  -  T°)/L  (Section  V) 
s  interface  position  along  the  diagonal,  m 
T  temperature,  K 

T  dimensionless  temperature,  (T0-T°)/(T^- T°)  (Section  V) 

T  dimensionless  temperature,  (T0-T°)/(Tjn-T°)  (Section  III) 
* 

T  transition  temperature,  K  (Section  IV) 

T  scaled  temperature,  T°-T°,  K  (Sections  III  and  V) 
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Nomenclature  (continued) 
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* 


k 


T°  temperature,  K 
T&  environment  temperature,  K 

Tc  reference  temperature  for  convection,  K  (Section  IV) 

Tc  reference  temperature  for  the  Clausius- Clapeyron  equation,  K  (Section 
VII) 

Tc  cold  surface  temperature,  °C  (Section  II) 

T£  dimensionless  cold  surface  temperature,  (T°  -  T°) / (T^  -  T°)  (Section 

V) 

T°  cold  surface  temperature,  K 

Tcr  reference  temperature  for  the  Clausius- Clapeyron  relationship,  K 
hot  surface  temperature,  °C  (Section  II) 

dimensionless  hot  surface  temperature,  (T^  -  T®)/(T^  -  T°)  (Section  V) 
hot  surface  temperature,  K 
initial  dimensionless  temperature,  (T?  -  T°) / (T^  -  T°) 

T°  initial  temperature,  K 

temperature  at  the  wick- wall  interface,  K 
T.  j  temperature  at  node  j-1  in  the  radial  direction  in  the  wick  region,  K 
Tj  temperature  at  last  node  near  liquid- vapor  interface  in  the  radial 
direction  in  the  wick  region,  K 
Tm  melting  temperature,  K 

T°  melting  (or  freezing)  temperature,  K 

Tq  reference  temperature,  K 

Tq  initial  temperature,  K  (Section  IV) 

TQ  reference  temperature,  °C  (Section  II) 

Tf  reference  temperature  for  radiation,  K 
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Nomenclature  (continued) 


T  saturation  temperature,  K  (Section  IV) 

T  saturation  temperature,  °C  (Section  II) 
s 

T®,T^  the  temperatures  defined  in  Fig.  5.2,  K 
Ty  outer  wall  surface  temperature,  K  (Section  VI) 

temperature  of  heating  surface  in  the  case  of  the  constant  heat 
flux,  °C  (Section  II) 

Ttf  dimensionless  wall  temperature,  (T°  -  T°)/(T^  -  T°)  (Section  V) 

T°  wall  temperature,  K 

t  time,  s 

U  axial  velocity  m/s  (Sections  IV  and  VII) 

UQ  inlet  velocity,  m/s 

U,V  dimensionless  velocities,  uH/a,  vH/a  (Section  II) 

U,V  dimensionless  velocities,  u/UQ,  v/UQ  (Section  III) 

U,V,V  dimensionless  velocities,  uH/a^,  vH/a^,  wH/a^  (Section  V) 
u,v  velocities,  m/s  (Sections  II  and  III) 
u,v,w  velocities,  m/s  (Section  V) 

V  net  velocity  vector,  m/s 

V  radial  velocity,  m/s  (Section  IV) 

VQ  radial  velocity  at  the  wall,  m/s 

v  radial  velocity,  m/s  (Section  VI) 

v^  radial  vapor  velocity  at  the  interface,  m/s 

*  * 
v  radial  velocity  based  on  the  guessed  pressure  p  ,  m/s 

V  width  of  the  thermal  cavities,  m  (Section  II) 

w  axial  velocity,  m/s  (Section  VI) 


Nomenclature  (continued) 


w  axial  velocity  based  on  the  guessed  pressure,  p  ,  m/s 
X,R  dimensionless  coordinate  directions,  x/D,  r/D  (Section  III) 

X,Y  dimensionless  coordinate  directions,  x/H,  y/H  (Section  II) 

X,Y,Z  dimensionless  coordinate  directions,  x/H,  y/H,  z/H  (Section  V) 
x,y,z  coordinate  directions  (Section  V) 
x,r  coordinate  directions  (Section  III) 
x,y  coordinate  directions  (Section  II) 

x  coordinate  in  the  axial  direction  (Sections  IV  and  VII) 

x,  •  initial  location  of  the  transition  region,  m 
t ,  1 

x  dimensionless  location  of  the  transition  region,  (x  -  x  .  )/ft 

t ,  l 

z  axial  cordinate,  m  (Section  VI) 

Greek  Symbols 

2 

a  thermal  diffusivity,  m  /s 

a  relaxation  factor  (Section  IV) 

P  coefficient  of  volumetric  thermal  expansion,  °C  1 
7  ratio  of  specific  heats,  Cp/cy 

T  fraction  of  transition  to  turbulent  flow 
6  wall  or  liquid- wick  thickness,  m 

£t  time  increment,  s 

6x  distance  between  nodes,  m 

261°  phase- change  temperature  range,  or  mushy  phase  range,  K  (Section 

III) 

2 £T°  phase-change  temperature  range,  (T^  -  1°),  K  (Section  V) 
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Nomenclature  (continued) 


ST*  £T°/(T°n-T°)  (Section  III) 

ST*  £T°/(t£-T°)  (Section  V) 

fl(T  -  T  )  Dirac  function 

Ar  radial  distance  between  nodes,  m 

AT  small  finite  temperature  interval  around  Tm  to  define  mushy 
e  porosity  (Sections  I  and  IV) 

e  emissivity  (Section  VI) 

$  dimensionless  temperature,  (T  -  T£)/(T^  -  T  )  or  (T  -  Tc)/T 

II) 

0  contact  angle  of  the  liquid,  deg  (Section  IV) 

«0  (T“  -  t°)/(t°  -  T°) 

X  length  of  mean  free  path,  m  (Section  IV) 

X  parameter  in  eqn.  5.20 

p  dynamic  viscosity,  kg/(m-s) 

v  kinematic  viscosity,  m  /s 

(  emissivity 

3 

p  density,  kg/m 

3 

p  reference  density,  kg/m 

3 

p ^  vapor  density  at  the  liquid- vapor  interface,  kg/m 

p'  density  correction 

* 

p  guessed  vapor  density 

* 

p  modified  density  (Section  I) 

a  Stefan- Boltzmann  constant,  V/(m^-K^) 

r  shear  stress  (Sections  TV  and  VII) 
r  dimensionless  time,  UQ  t/D  (Section  III) 


zone,  K 


i 

(Section 
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Nomenclature  (continued) 


2 

r  dimensionless  time,  a ^  t/H  (Section  V) 

<j>  general  variable  in  the  discretization  eqn.  (6.9) 
♦  viscous  dissipation  term  in  eqn.  (6.14) 
oi  surface  tension,  N/m  (Section  IV) 
w  porosity  (Section  VI) 

xr  =  3/4  'xr  =  1/4 


T 


Subscripts 
a  adiabatic 

B  "bottom"  neighbor  of  grid  P 

b  control- volume  face  between  P  and  B 

c  condenser 

cond  conduction 

conv  convection 

E  "east"  neighbor  of  grid  P 

e  control- volume  face  between  P  and  E  (Section  V) 

e  evaporator  or  control- volume  face  between  grid  P  and  E  (Section  VI) 
eff  effective 

f  transfer  fluid 

it  working  substance  in  the  liquid  state  in  the  wick 

fs  working  substance  in  the  solid  state  in  the  wick 

i  initial  condition,  or  inside  radius  of  the  pipe  (Section  III) 
i  initial  condition  (Section  V) 

i  vapor- liquid  interface  (Section  VI) 

in  inlet 
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Nomenclature  (continued) 


inf  melting  interface 

L  latent  heat 

l  liquid  PCM  or  latent  heat  (Section  III) 

i  liquid  where  there  is  liquid  motion  in  the  wick  (Section  IV) 
i  liquid  phase  (Section  V) 

l  liquid  or  liquid- wick  (Sections  I  and  VI) 

£e  wick  region  where  the  working  substance  is  in  the  liquid  state 
m  mass  or  mushy  phase  (Section  III) 
m  mushy  phase  (Section  V) 

m  wick  structure  material  (Section  VI) 

me  wick  region  where  the  working  substance  is  in  the  mushy  state 
N  "north"  neighbor  of  grid  p 

n  control- volume  face  between  P  and  N 

o  outer  surface  of  the  PCM  module  (Section  III) 
o  properties  at  the  liquid- vapor  interface  (Sections  I  and  IV) 

o  properties  of  the  injected  or  extracted  fluid  at  the  wall  (Section 

VII) 

P  grid  point 

p  PCM 

S  "south"  neighbor  of  grid  P 

s  solid  PCM  (Section  III) 

s  wick  structure  material  (Section  IV) 

s  control- volume  face  between  P  and  S,  or  solid  phase  (Section  V) 

s  saturation  (Section  VI) 

se  wick  region  where  the  working  substance  is  in  the  solid  state 
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Nomenclature  (continued) 


"top"  neighbor  of  grid  P 

control- volume  face  between  P  and  T 

vapor 

"west"  neighbor  of  grid  P 

control- volume  face  between  P  and  V  (Section 

wall 

derivative  with  respect  to  x 
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Section  I 


EXPERIMENTAL  AND  NUMERICAL  ANALYSIS  OF  LOV  TEMPERATURE 
HEAT  PIPES  WITH  MULTIPLE  HEAT  SOURCES 


1.1  SUMMARY 

A  numerical  analysis  and  experimental  verification  of  the  effects  of 
heat  load  distribution  on  the  vapor  temperature,  evaporation  and 
condensation  rates,  and  the  heat  transport  capacity  for  heat  pipes  with 
multiple  heat  sources  is  presented.  The  elliptic  conjugate  mass,  momentum 
and  energy  equations  in  conjunction  with  the  thermodynamic  equilibrium 
relations  and  appropriate  boundary  conditions  for  the  vapor  region,  wick 
structure,  and  the  heat  pipe  wall  are  given.  A  numerical  solution  of  these 
coupled  equations  as  a  conjugate  problem  was  required  for  the  performance 
prediction  of  heat  pipes  with  multiple  heat  sources.  The  experimental 
testing  of  a  copper- water  heat  pipe  with  multiple  heat  sources  was  also 
made  which  shows  excellent  agreement  with  the  numerical  results.  An 
optimization  of  the  heat  distribution  for  such  heat  pipes  was  performed  and 
it  was  concluded  that  by  redistribution  of  the  heat  load,  the  heat  capacity 
of  heat  pipes  can  be  increased. 
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1.2  INTRODUCTION 


Since  the  invention  of  the  heat  pipe  by  Grover  et  al.  (1964),  many 
investigations  have  been  performed  concerning  heat  pipe  operation,  factors 
limiting  heat  pipe  performance,  heat  pipe  applications  and  design 
modifications  for  the  improvement  of  heat  pipe  performance.  Because  of  the 
simplicity  of  design  and  ease  of  manufacture  and  maintenance,  heat  pipes 
have  found  applications  in  a  wide  variety  of  areas  (Chi,  1976;  Dunn  and 
Reay,  1982),  including  energy  conversion  systems,  cooling  of  nuclear  and 
isotope  reactors,  cooling  of  electronic  equipment  and  high  performance 
space  applications.  The  performance  of  a  heat  pipe  often  is  critical  and 
is  judged  in  part  by  the  amount  of  heat  a  unit  length  of  the  heat  pipe  can 
transport  under  a  uniform  heat  load. 

Advanced  spacecraft  electronic  cooling  systems  will  require  heat 
transport  loops  and  heat  pipes  which  have  multiple  heat  sources  in  series 
along  the  heat  flow  path  separated  by  uniform  and  nonuniform  distances.  In 
addition,  advanced  spacecraft  and  hypersonic  thermal  management 
requirements  also  motivate  the  analysis  and  performance  evaluation  of  heat 
pipes  with  multiple  heat  sources.  In  general,  the  analytical  and 
experimental  analyses  developed  for  single  evaporator  heat  pipes  cannot  be 
applied  to  multiple  evaporator  heat  pipes  with  nonuniform  heat  loads. 
Military  applications  of  such  heat  pipes  are  the  cooling  of  leading  edges 
and  nose  cones  of  re-entry  vehicles  and  cooling  rail  guns  and  laser 
mirrors.  Future  hypersonic  vehicle  structures  will  also  require  high 
performance  heat  pipes  with  multiple  evaporators  for  regions  subjected  to 
high  intensity  heating  where  large  quantities  of  heat  must  be  absorbed  and 
the  heating  distribution  is  not  uniform. 
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The  axial  heat  flux  in  heat  pipes  is  governed  in  principal  by  vapor 
flow  and  liquid  return  flow  limitations.  If  the  liquid  return  flow  can  be 
guaranteed  by  the  appropriate  wick  structure  design,  then  the  axial  heat 
flux  is  limited  only  by  the  vapor  flow  effects.  The  vapor  flow  is  limited 
by:  the  choking  phenomenon,  where  the  vapor  velocity  approaches  the  sonic 
velocity;  the  viscous  phenomenon,  in  which  the  vapor  pressure  approaches 
zero;  and  the  entrainment  phenomenon,  where  liquid  droplets  are  entrained 
into  the  vapor  flow. 

The  interruption  of  the  liquid  return  flow  can  be  caused  by 
insufficient  capillary  pressure  (wicking  limit)  or  by  the  formation  of 
bubbles  in  the  wick  (boiling  limit).  It  should  be  noted  that  the  capillary 
limitation  on  the  heat  transport  rate  (3  _  for  conventional  heat  pipes 
can  be  derived  from  the  heat  transport  factor  (QT)„  mo  if  the  heat  flow 
distribution  along  the  pipe  is  known.  In  conventional  heat  pipe  analysis 
(Chi,  1976),  the  heat  flux  distribution  is  assumed  uniform  for  simplicity 
and  because  there  is  no  information  available  for  nonuniform  heat  loading. 

Many  investigators  have  examined  the  problem  of  vapor  flow  in  circular 
and  annular  heat  pipes  under  uniform  heating  loads  using  numerical  and 
analytical  methods.  The  majority  of  these  investigators  neglected  the 
effect  of  the  wick  and  the  heat  conduction  in  the  pipe  wall.  Compressible 
vapor  flow  analysis  is  also  needed  for  the  sonic  conditions.  A  recent 
detailed  literature  survey  concerning  the  thermal  modelling  of  heat  pipes 
is  given  by  Chen  and  Faghri  (1989).  The  review  covers  both  analytical  and 
numerical  methods,  transient  and  steady  conditions  as  well  as  one-  or 
two-dimensional  modeling.  Some  of  the  most  recent  investigations  are  by 
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Busse  (1987),  Faghri  and  Parvani  (1988),  Issacci  et  al.  (1988),  Bowman  and 
Hitchcock  (1988),  Cao  et  al.  (1989b),  Faghri  and  Thomas  (1989),  Faghri 
(1989) ,  Chen  and  Faghri  (1989),  Faghri  and  Chen  (1989),  Jang  et  al. 

(1989a),  and  Jang  et  al.  (1989b). 

Gernert  (1986)  attempted  the  analysis  of  heat  pipes  with  multiple  heat 
sources  with  the  use  of  superposition  and  the  extension  of  the  existing 
theories  for  a  single  evaporator  heat  pipe.  It  is  clear  that  one  cannc. 
use  superposition  for  fluid  flow  analysis  due  to  the  nonlinearity  of  the  « 

momentum  equations.  Furthermore,  axial  conduction  is  neglected  which  can 
play  an  important  role  with  multiple  heat  sources.  The  objective  of  this 
paper  was  to  achieve  the  following: 

1-  The  formulation  of  a  detailed  mathematical  model  with  numerical 
results  for  predicting  the  performance  of  heat  pipes  with  multiple 
evaporators.  This  was  accomplished  by  solving  numerically  the 
complete  mass,  momentum  and  energy  equations  in  the  vapor  region,  wick 
structure  and  the  heat  pipe  wall  for  heat  pipes  with  multiple  heat 
sources  as  a  conjugate  problem.  This  new  model  includes  the  effect  of 
liquid  flow  in  the  wick  which  is  important  for  low  temperature  heat  * 

pipes  rather  than  assuming  pure  conduction  through  the  wick  as  was 
done  by  Chen  and  Faghri  (1989).  t 

2  The  experimental  testing  of  a  copper- water  heat  pipe  with  multiple 
evaporators  using  various  heat  loads  in  each  evaporator  for  model 
verification  as  well  as  optimization  of  the  heat  distribution  for  such 
heat  pipes. 
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1.3  MATHEMATICAL  MODELING  AND  NUMERICAL  ANALYSIS 


Fluid  mechanics  and  heat  transfer  problems  in  heat  pipes  are 
categorized  in  the  four  basic  types:  vapor  flow  in  the  core  region;  liquid 
flow  in  the  wick;  heat  conduction  in  the  heat  pipe  wall;  and  the 
interaction  between  the  liquid  and  vapor  flows.  Most  of  the  analytical  and 
numerical  studies  on  heat  pipes  have  been  done  on  the  vapor  region  for  the 
case  of  uniform  heating  and  cooling  with  one  evaporator  section. 

The  physical  problem  under  consideration  is  a  conventional  heat  pipe 
with  multiple  heat  sources  as  shown  in  Fig.  1.1a,  which  is  divided  into 
three  regions,  namely,  the  vapor  flow  region,  the  liquid- wick,  and  the 
wall.  It  is  assumed  that  the  vapor  flow  in  all  segments  of  the  heat 
pipe,  i.e.,  the  evaporator,  adiabatic  and  the  condenser  sections,  is 
operating  under  laminar,  subsonic  and  steady  conditions.  Nonuniform  radial 
inflow  and  outflow  boundary  conditions  are  needed  to  model  evaporation  and 
condensation.  Since  the  heat  pipe  is  closed  at  both  ends,  it  is  required 
that  the  vapor  which  flows  out  of  the  evaporator  segment  enters  into  the 
condenser  section.  The  conservation  of  mass,  momentum,  and  energy 
equations  for  the  compressible  vapor  flow  region  including  viscous 
dissipation  are: 


£-(p  w  )  +  -  p  rv  )  =  0 
ozxvv  v;  r  5rVPv  v' 


(1.1) 
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Figure  1.1a.  The  Multiple  Evaporator  Heat  Pipe  and  Coordinate  System 
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•  Wall  Thermocouple 


Figure  1.1b.  Low  Temperature  Heat  Pipe  with  Multiple  Heat 

Source  Thermocouple  locations 
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It  should  be  noted  that  the  perfect  gas  law  is  employed  to  account  for  the 
compressibility  of  the  vapor. 


The  use  of  liquid  capillary  action  is  the  unique  feature  of  the  heat 
pipe.  From  a  fundamental  viewpoint,  the  liquid  capillary  flow  in  heat 
pipes  with  screen  wicks  should  be  modeled  as  a  flow  through  a  porous  media. 
It  is  assumed  that  wicking  material  is  isotropic  and  of  constant  thickness. 
In  addition,  the  wick  is  saturated  and  the  vapor  condenses  or  the  liquid 
evaporates  at  the  vapor- liquid  interface. 

The  averaging  technique  was  applied  by  many  investigators  to  obtain 
the  general  governing  equation  which  describes  the  conservation  of  momentum 
in  a  porous  structure.  Since  the  development  of  Darcy’s  semi- empirical 
equation,  which  characterizes  the  fluid  motion  under  certain  conditions, 
many  scientists  have  tried  to  develop  and  extend  Darcy’s  law  in  order  to 
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see  the  effect  of  the  inertia  terms.  In  this  respect,  those  who  have  tried 
to  model  the  flow  with  the  Navier- Stokes  equation  were  the  most  successful. 


The  general  equations  of  motion  and  energy  for  steady  state  laminar 
incompressible  liquid  flow  in  porous  media  in  terms  of  the  volume- averaged 
velocities  (Bachmat  and  Bear,  1986)  are: 


dvi  dvi  * i  rid, 

pivi  w  +  pivi  w  -  ■  w  +  pe  [  t  ie  (r  +  ~^r  J  — C1-5) 


dv,  dv  a  dp.  r  a  i  5(v/r)  d2yf  i 

pivi  W  +  pivi  W  =  '  W  +  H  [  3F  (r  “3r — >  +  ^2~  . 


We  (A  ,, 
— K —  (l-6) 


81.  81.  ri  B  81.  d2 T,  , 

^Cp,f  (vf  <?r""  +  wf  cTz-)  keff  [r  3F  (r  +  j 


(1.7) 


For  only  a  capillary  jacket  when  the  pores  in  the  isotropic  wicking 
material  are  completely  filled  with  liquid, 


c 


r 


K 


r 


=  K 


(1.8) 


The  above  assumption  was  made  in  all  the  numerical  results  presented  in 
this  paper.  The  method  of  calculating  the  effective  thermal  conductivity 
of  the  wick,  k  from  the  thermal  conductivity  of  the  solid  and  liquid 
phases  are  given  by  Chi  (1976), 
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The  steady  state  energy  equation  that  describes  the  temperature  in  the  heat 
pipe  wall  is: 


d 


ft  i  a  f  ft.  ■ 
kw  w\  +  r  tTr  (kwr3F~. 


=  0 


(1.10) 


where  ky  is  the  local  thermal  conductivity  of  the  heat  pipe  wall.  The 
boundary  conditions  are: 

w(0,r)  =  v(0,r)  =  0 

w(LT,r)  =  v(LT,r)  =  0 

w(z,a)  =  0 


w(z,b)  =  0 


(l.n) 

(1.12) 

(1.13) 

(1.14) 
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b  * 

"keff  W 
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v  fg 


v(z,b)  =  0 
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T(z,a)  = 
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(1.18) 


In 


h 

p0 


dT  n 

<IF  (z’c)  =  I  (z)  (1 . 19) 

In  boundary  condition  (1.19),  the  heat  flux  was  constant  under  the 
active  evaporators  and  zero  in  all  adiabatic  and  transport  sections  and  the 
inactive  evaporators.  In  the  condenser  section,  the  heat  flux  was  assumed 
to  be  uniform  based  on  the  total  heat  input  through  the  evaporators.  The 
temperature  along  the  vapor- liquid  interface  is  taken  as  the  equilibrium 
saturation  temperature  corresponding  to  its  equilibrium  pressure  condition. 
By  the  above  procedure,  the  effect  of  energy  conduction  in  the  heat  pipe 
wall  and  the  liquid  flow  in  the  wick  is  taken  into  account. 

For  simplicity  and  generality,  the  problem  should  be  solved  using 
conjugate  heat  transfer  analysis  as  a  single  domain  problem.  To  achieve 
this,  one  should  generalize  the  conservation  equations  for  mass,  momentum 
and  energy  such  that  each  conservation  equation  should  have  the  same  source 
term  for  the  energy  equation  in  terras  of  temperature  as  a  dependent 
variable  for  the  three  regions  (i.e.,  wall,  wick  and  vapor).  In  addition, 
the  continuity  of  temperature,  heat  flux  and  mass  flux  should  be  satisfied 
at  each  of  the  interfaces  as  well  as  some  special  boundary  conditions  like 
the  no- slip  condition,  thermodynamic  equilibrium,  etc.  The  following 
approach  is  not  only  very  beneficial  for  those  who  develop  their  own 
computer  program  to  solve  this  heat  pipe  problem,  but  also  for  those  who 
wish  to  use  the  existing  commercial  codes  to  set  up  this  problem. 
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The  energy  equation  can  be  written  in  terms  of  temperature  as 
DT  4.  S  /<  n/-.\ 

p  ut  -  r  v  T  +  r  t1-20) 

p  p 

Here  S  is  the  source  term  which  includes  viscous  dissipation  and  pressure 
work.  It  should  be  noted  that  solving  the  problem  in  terms  of  enthalpy 

( 

does  not  preserve  the  condition  of  temperature  continuity  at  an  interface 
when  harmonic  averaging  is  employed  and  leads  to  an  incorrect  calculation 

% 

of  diffusion. 


For  the  vapor  region,  the  energy  equation  is 

(1. 

Multiplying  both  sides  of  the  energy  equation  for  the  liquid- wick  region  by 

c  Jc  results  in 
p,r  p,v 


DT 
\  Dt 


—  ?2t  +  A 


P,v 


P,v 


n  ^  K  =  p-  '2t  +  r 

p,v  p,v 


P,V 


(1.22) 


Similarly,  the  energy  equation  for  the  heat  pipe  wall  is 


Cp,w  DT  ,.,N  y2T  _S 

"  CP,*  m  V  Cp,T 


(1.23) 


From  this  transformation,  we  observe  the  following: 
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1)  The  source  term  for  the  energy  equation  is  dividea  by  Cp  y  for  all 
three  regions. 

* 

2)  One  needs  to  set  up  the  density  equal  to  the  modified  density  p  for 
different  regions,  i.e.,  p  for  vapor,  pf  c  Jc  for  liquid  and  p 

cp,»/cp,v  for  solid- 

3)  The  diffusion  coefficients  contain  only  one  c  ,  namely,  c  .  The 

P  P  ■> v 

transformation  makes  possible  an  exact  representation  of  diffusion 
across  an  interface  when  the  harmonic  average  is  used  for  the 
diffusion  coefficient  at  the  interfaces. 

The  momentum  equation  for  the  vapor  flow  does  not  need  any  special 
* 

treatment  since  p  =  p^  in  the  vapor  region.  In  the  wick  region,  the 

momentum  equation  can  be  transformed  in  the  same  manner. 


*  DV 
p  Bt  - 


* 

Vp  +  vt  p  V  - 


We  observe  the  following: 


(1.24) 


1)  A  source  term  should  be  added  in  the  momentum  equation  for  the 
porosity  effect  in  the  wick  region. 

* 

2)  The  pressure  solved  for  is  the  modified  pressure  p  ,  which  is 

proportional  to  the  actual  pressure.  The  actual  pressure  drop  between 

* 

two  points  in  the  flow  fluid  can  be  calculated  as  Ap  =  c  /c  »  Ap 

P  J  V  P  }  C 

The  transformed  continuity  equation  can  be  written  as 
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(1.25) 


V  •  (p*V)  =  0 

The  finite- difference  iteration  method  of  solution  developed  by 
Spalding  (1980)  is  employed  in  the  solution  of  the  elliptic  governing 
equations  (1.1-1.10)  subjected  to  the  boundary  conditions  (1.11-1.19).  In 
this  methodology,  finite- domain  equations  are  derived  by  integration  of  the 
differential  equations  over  a  control  volume  surrounding  a  grid  node.  The 
source  terms  due  to  viscous  dissipation  and  pressure  work  in  the  energy 
equation,  and  the  source  term  due  to  the  porous  matrix  in  the  momentum 
equation  are  linearized,  and  the  "SIMPLEST"  practice  (Spalding,  1980)  was 
employed  for  the  momentum  equations.  The  solution  procedure  is  based  on  a 
line-by-line  iteration  method  in  the  axial  direction  and  the  Jacobi 
point- by- point  procedure  in  the  radial  direction. 


The  energy  equation  is  not  continuous  at  the  vapor- liquid  interface 


due  to  the  latent  heat  of  evaporation  and  condensation.  The  term  m  h^ 
must  be  added  as  a  source  term  in  the  energy  equation  at  the  vapor- liquid 
interface.  For  the  first  few  iterations,  it  was  assumed  the  heat  flux  at 
the  liquid- vapor  interface  is  equal  to  the  outer  wall  heat  flux.  An  exact 
energy  balance  given  by  eqn.  (1.15)  is  satisfied  after  these  initial 
iterations. 


The  computation  proceeded  until  the  sum  of  the  absolute  volumetric 
errors  over  the  whole  field  was  negligibly  small  (<10  ).  The  numerical 
results  were  also  tested  for  grid  independence  by  systematically  varying 
the  number  of  grids  in  both  the  r-  and  z- directions.  For  the  numerical 
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results  presented  here,  the  final  grid  sizes  for  all  of  the  cases  presented 
were  chosen  as  follows:  100  x  (20  +  5  +  10)  =  (axial)  x  [(radial  vapor)  + 

(radial  liquid-wick)  +  (radial  solid  wall)]. 

1.4  EXPERIMENTAL  APPARATUS  AND  PROCEDURE 

The  objective  of  this  project  was  to  determine  the  operating 

characteristics  and  performance  of  a  low  temperature  heat  pipe  with 
multiple  heat  sources.  The  low  temperature  heat  pipe  was  a  copper- water 
heat  pipe  designed  to  operate  at  a  vapor  temperature  of  60-100°C.  The  * 

compatibility  and  efficiency  of  the  copper- water  heat  pipe  is 

well- documented  in  heat  pipe  literature. 

The  low  temperature  heat  pipe  shell  and  end  caps  were  fabricated  from 
standard  composition  copper  (UNS- C12200) .  The  heat  pipe  shell  is  1000  mm 
in  length  with  an  outside  diameter  of  25.4  mm,  and  a  wall  thickness  of 
1.7  mm.  The  ends  of  the  heat  pipe  shell  were  machined  so  that  the  end  caps 
fit  snugly.  A  simple  circumferential  screen  wick  consisting  of  two  wraps 
of  50  per  inch  mesh  copper  screen  was  installed  to  provide  a  liquid  return 
path  to  the  evaporators.  A  summary  of  the  design  specifications  for  the 
heat  pipe  is  given  in  Table  1.1.  * 

All  of  the  heat  pipe  parts  were  carefully  fitted,  cleaned,  and  ♦ 

deoxidized  using  standard  procedures  (Chi,  1976),  and  the  end  caps  were 

brazed  to  the  heat  pipe  shell  using  Harris  15  filler  rod.  The  interior  of 

the  heat  pipe  was  protected  from  oxidation  during  brazing  with  a  cover 
layer  of  hydrogen.  A  bellows  type  valve  was  attached  to  the  heat  pipe  fill 
tube  to  facilitate  sealing,  purging,  and  charging  of  the  fluid  inventory 
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TABLE  1.1 


DESIGN  SUMMARY  OF  THE  LOV  TEMPERATURE  HEAT  PIPE 
VITH  MULTIPLE  HEAT  SOURCES 


Container  material 

copper 

Vick  material 

copper 

End  cap  material 

copper 

Length 

1000  mm 

Container  O.D. 

25.4  mm 

Container  I.D. 

22.0  mm 

End  cap  thickness 

3.175  mm 

Screen  mesh  number 

1.97  x  103  m  1 
(50  mesh) 

Screen  wire  diameter 

.178  mm 

Screen  wick  thickness 

.712  mm 

Vapor  core  diameter 

20.5  mm 

Working  fluid 

distilled  water 

Fluid  charge 

40  cm3 

* 

Evaporator  length 

(4)  @  63.5  mm 

Transport  length 

180  mm 

** 

Condenser  length 

300  mm 

Evaporators  separated  by  75.3  mm  adiabatic  sections,  with  a  20  mm 
adiabatic  section  between  the  evaporator  end  cap  and  evaporator  1. 

20  mm  adiabatic  section  between  the  condenser  and  the  condenser  end 
cap. 
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inside  the  heat  pipe. 


The  heat  pipe  has  four  independently  controlled  evaporators  in  the 
evaporator  section.  Each  evaporator  consists  of  a  63.5  mm  long  257  watt 
Minco  Thermofoil  heater  clamped  to  the  heat  pipe  (Fig.  1.1b).  Power  input 
for  each  heater  was  supplied  by  a  120  volt  variable  ac  transformer  and 
measured  with  a  Fluke  77  multimeter  which  has  an  accuracy  of  approximately 
±27.  of  the  reading.  Heat  was  removed  from  the  condenser  section  by  a  300 
mm  long  copper  water-cooled  manifold  calorimeter  mounted  to  the  heat  pipe. 
Thermocouples  in  the  calorimeter  inlet  and  outlet  and  mass  flow 
measurements  were  used  to  calculate  the  power  output  from  the  calorimeter, 
which  was  compared  to  the  electrical  input  power  so  that  accurate  heat  pipe 
power  levels  were  obtained. 

All  heat  pipe  vapor  and  wall  temperatures  were  measured  with  standard 
Type  K  thermocouples,  which  have  an  accuracy  of  ±2.2°C  or  0.757,  of  the 
reading,  whichever  is  greater.  There  were  thirteen  wall  thermocouples  and 
six  vapor  space  thermocouples  (Fig.  1.1b).  Vail  temperatures  in  the 
adiabatic  sections  were  measured  with  beaded  thermocouples  clamped  to  the 
heat  pipe  shell.  Four  sheathed  thermocouples  were  soft- soldered  to  the 
heat  pipe  wall  within  the  calorimeter.  These  thermocouples  exit  the 
calorimeter  through  shallow  grooves  machined  in  the  heat  pipe  wall.  Vail 
temperatures  at  the  axial  center  of  each  evaporator  were  measured  with 
sheathed  thermocouples  soft- soldered  to  the  heat  pipe  shell.  These 
thermocouples  exit  the  evaporator  sections  through  a  shallow  groove  in  the 

_  7 

heat  pipe  wall.  The  heat  pipe  was  evacuated  to  a  pressure  of  10  torr 
prior  to  wall  thermocouple  soldering  to  prevent  interior  oxidation.  Vapor 
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temperatures  were  measured  with  a  4.7  mm  diameter  type  316  stainless  steel 
sheathed  multipoint  thermocouple  mounted  axially  within  the  heat  pipe  vapor 
space.  The  multipoint  thermocouple  exits  the  heat  pipe  through  a  tube  in 
the  evaporator  end  cap  and  was  swaged  in  place  to  provide  a  leak- tight 
seal. 


'  After  the  installation  of  all  thermocouples  and  heaters,  the  heat  pipe 

was  processed  in  a  specially  built  heat  pipe  filling  station.  The  heat 
t  pipe  was  evacuated  to  a  pressure  of  10  torr  and  filled  with  40  cm  of 

degassed,  distilled  water. 

A  schematic  of  the  test  setup  is  shown  in  Fig.  1.2.  The  heat  pipe  was 
mounted  horizontally  on  an  optical  bench.  The  optical  bench  was  equipped 
with  adjustable  base  plates  to  allow  for  precise  leveling  of  the  heat  pipe 
assembly.  The  thermocouples  were  read  by  a  Fluke  2285B  Data  Logger. 
Cooling  water  was  supplied  to  the  calorimeter  by  a  centrifugal  pump 
connected  to  a  constant  head  reservoir  in  order  to  maintain  a  steady 
coolant  flow.  A  coolant  pre- heater  was  located  prior  to  the  calorimeter 
inlet  to  allow  greater  flexibility  in  the  heat  pipe  operating  temperature. 

*  The  entire  length  of  the  heat  pipe  assembly  was  insulated  with  2  inches  of 
Fiberfrax  ceramic  fiber  insulation. 

* 

Tests  were  performed  to  establish  the  maximum  heat  transport  capacity 
of  the  heat  pipe  with  variable  heat  input  rates  and  locations  when  the  heat 
pipe  was  horizontal.  The  operating  vapor  temperature  was  held  between  65°C 
and  85°C.  Energy  balances  between  the  heat  input  by  the  heaters  and  the 
heat  removed  by  the  calorimeter  were  monitored  to  ensure  an  energy  balance 
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of  at  least  907.  during  steady  state  operation. 


The  primary  interest  in  these  tests  was  to  determine  the  capillary 
limit  of  the  heat  pipe  in  multiple  evaporator  operation.  Because  the 
evaporator  located  farthest  from  the  condenser  will  experience  dry- out 
first,  all  multiple  evaporator  tests  were  performed  with  some  combination 
of  evaporator  1  and  the  other  evaporators.  The  maximum  heat  transport 
capacity,  or  capillary  limit,  of  the  heat  pipe  was  defined  as  a  sudden, 
rapid,  and  continuous  increase  in  the  wall  temperature  of  evaporator  1. 

1.5  RESULTS  AND  DISCUSSION 

The  heat  pipe  was  first  operated  in  the  single  evaporator  mode  in 
order  to  establish  the  capillary  limits  of  the  heat  pipe  at  all  four  heat 
input  locations.  The  capillary  limit  for  the  evaporator  farthest  from  the 
condenser,  evaporator  1,  was  found  to  be  181  watts.  Limits  for  evaporators 
2,  3,  and  4  were  greater  than  the  maximum  power  levels  that  the  heaters 
could  deliver,  which  was  250  watts.  Figure  1.3  shows  the  difference 
between  the  evaporator  wall  temperature  and  the  vapor  temperature 
(ATwaH_ vapor)  versus  t^ie  Power  input  for  each  single  evaporator  test.  The 
single  evaporator  tests  show  that  the  ATwa^  vap0r  increases  linearly  with 
increasing  power,  and  that  the  heat  pipe  can  successfully  operate  with  a 
high  AT 

°  wall- vapor 

The  test  matrix,  Table  1.2,  shows  the  experimentally  determined 
capillary  limits  for  the  heat  pipe  in  the  single  and  multiple  evaporator 
operating  modes.  All  of  the  multiple  evaporator  tests  were  performed  with 
evaporator  1  in  combination  with  the  other  evaporators,  since  the  single 
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Figure  1.3.  ^Twa||— vap0r  versus  heat  input  for  different  single  evaporator  operations 


TABLE  1.2 


LOV  TEMPERATURE  BEAT  PIPE  WITH  MULTIPLE 
HEAT  SOURCES  TEST  MATRIX 


Power  listed  is  experimental  capillary  limit++ 


Power  in  Evaporator  (watts) 


Case 

1 

2 

3 

4 

Total  Power 

1 

181 

- 

- 

- 

181 

2 

- 

>247  * 

- 

- 

>247  * 

3 

- 

- 

>246  * 

- 

>246  * 

4 

- 

- 

- 

>246  * 

>246  * 

5 

119 

118 

- 

- 

237 

6 

112 

112 

110 

- 

334 

7 

97 

98 

99 

98 

392 

8 

161 

- 

160 

- 

321 

9 

150 

- 

- 

149 

299 

10 

97  ** 

- 

- 

>245  * 

>342  * 

11 

100  ** 

- 

250 

- 

350 

12 

102  ** 

182 

_ 

_ 

284 

The  capillary  limit  for  multiple  evaporator  operation  was  defined  as  a 
sudden,  rapid  and  continuous  increase  in  the  wall  temperature  in 
evaporator  1. 

* 

The  capillary  limit  was  not  reached  in  cases  2,  3,  4,  and  10  due  to 
the  maximum  power  level  limitations  of  the  heaters. 

In  cases  10,  11,  and  12  evaporator  1  was  held  constant  at  100  watts 

(nominal) . 
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evaporator  tests  indicated  that  the  evaporator  located  farthest  from  the 
condenser  had  the  lowest  capillary  limit.  Figures  1.4  -  1.7  show  typical 
plots  of  wall  and  vapor  temperatures  versus  axial  location  for  the  heat 
pipe  in  multiple  evaporator  operation  for  cases  1,  5,  6,  and  7, 
respectively,  at  power  loads  significantly  less  than  the  capillary  limit. 
The  numerical  results  based  on  the  mathematical  modeling  presented  in  this 
paper  are  also  shown  in  Figures  1.4- 1.7  which  show  good  agreement  with 
experimental  data  for  the  wall  and  vapor  temperatures  along  the  heat  pipe. 
In  the  condenser  section,  it  can  be  seen  that  the  experimental  and 
numerical  values  of  the  outer  wall  temperature  do  not  coincide.  This  is 
due  to  the  assumption  of  a  constant  heat  flux  along  the  condenser  section 
at  the  outer  wall,  which  is  not  precisely  correct  in  this  particular 
design,  and  the  fact  that  the  thermocouples  mounted  on  the  condenser  wall 
may  be  affected  by  the  cooling  water  circulating  in  the  calorimeter.  .  Near 
the  evaporator  end  cap,  the  outer  wall  temperature  can  be  seen  to  decrease 
from  the  center  of  evaporator  1  to  the  end  cap.  The  reason  for  this  is  the 
20  mm  unheated  length  between  the  evaporator  end  cap  and  evaporator  1.  It 
should  be  noted,  however,  that  even  though  the  outer  pipe  wall  is  adiabatic 
in  the  unheated  lengths,  these  sections  still  act  as  evaporators  due  to 
axial  conduction. 

In  Figs.  4. 5-4. 7,  the  experimental  wall  temperatures  in  the  center  of 
the  active  evaporators  decrease  closer  to  the  condenser,  but  the  numerical 
model  predicts  that  these  peak  temperatures  are  all  nearly  the  same. 
Again,  this  is  due  to  the  assumption  of  a  constant  heat  flux  at  the  outer 
wall  in  the  condenser  section.  Since  the  heat  pipe  is  a  closed  system,  the 
downstream  boundary  conditions  can  affect  the  temperature,  velocity,  and 
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■  Experimental  Vapor  Temp.  Numerical  Wall  Temp. 

□  Experimental  Wall  Temp.  Numericd  Vapor  Temp. 


Figure  1.6.  Heat  pipe  wall  and  vapor  temperatures  versus  axial  location  for  Case  6 


density  upstream. 


In  cases  5-9  the  heat  pipe  was  operated  with  uniformly  increasing 
power  to  each  evaporator  until  dry- out  occurred  in  evaporator  1.  Cases  5, 
6,  and  7  indicate  that  the  total  heat  load  on  the  heat  pipe  can  be 
significantly  increased  by  adding  evaporators,  although  the  maximum  heat 
flux  for  evaporator  1  decreases  as  evaporators  are  added  downstream  (i.e., 
toward  the  condenser  section) . 

In  cases  10- 12  the  power  input  to  evaporator  1  was  held  at  a  constant 
100  watts,  while  the  power  input  to  evaporators  2,  3  or  4  was  increased 
uniformly.  Figures  1.8  and  1.9  show  that  vap0r  evaporator  1 
remains  relatively  constant  as  the  total  heat  load  on  the  heat  pipe  is 
increased  for  cases  11  and  12,  respectively,  then  evaporator  1  suddenly 
drys  out  when  a  specific  downstream  power  level  is  reached.  This  maximum 
downstream  power  level  decreases  as  the  second  evaporator  moves  farther 
away  from  the  condenser. 

Bienert  et  al.  (1977)  demonstrated  a  jet  pump- assisted  arterial  heat 
pipe  in  which  a  venturi  located  in  the  vapor  space  maximized  capillary 
pumping  by  providing  a  low  pressure  source  to  the  evaporator  end  of  the 
heat  pipe  artery.  Similarly,  the  addition  of  one  or  more  evaporators  near 
the  condenser  of  a  heat  pipe  will  provide  a  pressure  drop  in  the  wick  and 
increase  the  liquid  mass  flow  rate,  while  the  loss  of  capillary  pressure  in 
the  wick  due  to  the  additional  evaporators  is  minimal  given  their  proximity 
to  the  condenser.  This  argument  can  be  applied  to  cases  8  and  9  that  show 
the  addition  of  one  downstream  evaporator  increases  the  heat  capacity  of 
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the  heat  pipe  by  407.  over  the  single  evaporator  case  1.  It  is  also  true 
that  in  the  case  of  multiple  heat  sources,  one  is  reducing  the  effective 
length  of  the  heat  pipe  and  therefore  increasing  the  total  heat  capacity. 


1.6  CONCLUSIONS 

A  complete  two-dimensional  numerical  modelling  of  heat  pipes  with 
multiple  heat  sources  including  the  heat  conduction  in  the  wall,  fluid  flow 
is  a  porous  media  for  the  screen  wick  as  well  as  compressible  vapor  flow  as 
a  conjugate  problem  is  given.  This  methodology  can  be  easily  extended  to 
predict  the  performance  of  capillary  pumped  loops  as  well. 

A  simple  circumferential  screen- wick  copper- water  heat  pipe  with 
multiple  heat  sources  has  also  been  fabricated  and  successfully  tested. 
These  tests  show  that  the  maximum  total  heat  load  on  the  heat  pipe  varies 
greatly  with  the  location  of  the  local  heat  fluxes.  Significant  heat  loads 
can  be  carried  by  evaporators  close  to  the  condenser  without  affecting  the 
operation  of  upstream  evaporators,  but  adding  additional  evaporators  to  the 
heat  pipe  lowers  the  maximum  heat  flux  capacity  of  the  evaporator  located 
farthest  from  the  condenser.  The  evaporator  wall- vapor  temperature 
difference  seems  to  be  mainly  a  function  of  the  local  heat  flux,  and  cannot 
be  used  to  predict  the  maximum  heat  capacity  of  the  heat  pipe  operating  in 
the  multiple  evaporator  mode.  The  numerical  results  agree  with  the 
experimental  results  for  wall  and  vapor  temperatures  along  the  heat  pipe 
for  cases  with  different  heat  distributions  along  the  evaporator  section. 
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Section  II 


HEAT  TRANSFER  IN  LIQUID  METALS  BY  NATURAL  CONVECTION 


2.1  SUMMARY 

A  numerical  study  of  the  natural  convection  in  low  Prandtl  number 
fluids  which  are  heated  from  below  has  been  performed.  Because  of  the 
limitations  of  experiments,  a  numerical  study  is  needed  to  obtain 
temperature  distributions  prior  to  boiling  and  to  study  mechanisms  of 
boiling  incipience  on  the  wall  for  fluids  with  very  low  Prandtl  numbers. 
The  finite- difference  scheme  SIMPLER  was  used  to  solve  the  steady 
two-dimensional  governing  equations  with  the  Boussinesq  approximation.  The 
numerical  results  were  compared  to  previous  experimental  results  with  an 
excellent  agreement. 


\ 
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2.2  INTRODUCTION 


Natural  convection  in  horizontal  enclosures  with  heating  from  below 
has  been  the  subject  of  many  numerical  studies  in  recent  years.  However, 
little  attention  has  been  paid  to  the  fluids  having  very  low  Prandtl 
numbers  such  as  liquid  metals  (McDonough  and  Catton,  1982;  Raithby  and 
Hollands,  1985;  Yang,  1987;  Reed,  1987).  McDonough  and  Catton  (1982) 
studied  natural  convection  using  a  mixed  finite  difference  -  Galerkin 
procedure  in  a  horizontal  square  box  which  was  heated  from  below,  cooled  * 

from  above  and  the  side  walls  were  insulated.  They  found  that  the 
numerical  results  with  lower  Pr  did  not  converge  as  quickly  as  those  with  * 

higher  Pr,  and  the  convergence  is  not  monotonic.  They  believed  that  this 
is  due  to  the  increasing  nonlinearity  as  the  Prandtl  number  decreases.  In 
the  numerical  results  presented  by  McDonough  and  Catton  (1982) ,  the  Prandtl 
numbers  were  not  lower  than  Pr  =  0.71.  Obviously,  further  numerical 
studies  in  the  low  Prandtl  number  range  are  needed,  which  is  one  of  the 
objectives  of  this  paper. 

For  pool  boiling  and  two-phase  flow  heat  transfer,  the  temperature 
distribution  in  liquid  metals  is  crucial  to  the  operation  and  boiling 
incipience.  In  a  gravitational  environment,  heat  transfer  in  liquid  metals 
prior  to  boiling  incipience  is  a  problem  of  natural  convection  combined  4 

with  conduction,  even  for  liquid  metal  layers  with  a  thickness  on  the  order 
of  10  mm.  Because  of  the  limitations  of  experiments,  a  numerical  study  is  i 

needed  to  obtain  temperature  distributions  prior  to  boiling  and  to  study 
mechanisms  of  boiling  incipience  on  the  wall  for  very  low  Prandtl  numbers. 

This  is  another  motivation  of  the  present  numerical  study. 
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Many  pool  boiling  and  evaporation  test  sections  for  liquid  metal 
layers  can  be  modeled  as  a  two-dimensional  partial  heating  configuration  as 
shown  in  Fig.  2.1,  with  the  dimension  perpendicular  to  the  paper  being 
infinitely  long.  A  uniform  heat  flux  along  the  heating  element  is 
specified.  The  upper  liquid  surface  is  kept  at  the  saturation  temperature 
Tg  =  Tc,  corresponding  to  the  system  pressure,  and  is  considered  as  a  free 
surface  in  contrast  with  the  rigid  lower  surface.  The  container  wall 
except  the  heating  element  is  insulated  to  prevent  heat  loss  to  the 
environment.  Another  alternative  boundary  condition  for  this  problem  is  to 
specify  a  constant  temperature  T^  >  Tg  at  the  lower  surface,  which  has  been 
used  more  often  in  the  previous  numerical  studies. 

In  this  paper,  the  natural  convection  of  fluids  having  very  low 
Prandtl  numbers  down  to  Pr  =  0.00125  in  enclosures  with  partial  or  full 
heating  from  below  will  be  studied  numerically,  and  the  numerical  results 
will  be  compared  with  the  existing  experimental  data  (Takenaka  et  al., 
1983)  and  the  empirical  heat  transfer  equation  (Raithby  and  Hollands, 
1985),  respectively. 

2.3  FORMULATION  AND  NUMERICAL  METHOD 

The  dimensionless  governing  equations  for  steady  two-dimensional 
laminar  flow  of  a  Newtonian  fluid  with  no  dissipation  under  Boussinesq 
approximation  can  be  written  as  follows: 


Fig.  2.1  Two-dimensional  horizontal  liquid  metal 
layer  with  localized  heating  from  below 
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with  the  following  dimensionless  variables: 


x  y  H  H 

X  =  H’Y  =  H’U  =  ua;V  =  va;^  =  T  ’ 


T  -  T 


P  =  (p  +  P0 gy)  H 2 /pa2  ;  Pr  =  v/a  ;  Raj  =  g/ffl3T Jav 


where 


"o )/<T  -  To> 

For  the  partial  heating  from  below,  two  geometrical  parameters  are 
needed: 


Bj  =  V/H,  and  B2  =  L/H. 

The  boundary  conditions  for  the  constant  heat  flux  and  free  upper 
surface  shown  in  Fig.  2.1  are 
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5U 

V  =  0  »  =  0  ,  0  =  0  Y  =  1  ,  -B1<X_<B1  (2.5a) 

80 

V  =  U  =  0,^  =  0  0  <  Y  <  1  ,  X  =  -Bj  (2.5b) 

80 

V=U=0,^=0  0  <  Y  <  1  ,  X  =  Bj  (2.5c) 

For  the  rigid  lower  boundary  condition  (Y  =  0)  t 

V  =  U  =  0  -B1  <  X  <  B1  (2.5d)  t 

80  q  H 

gy  =  ~  'jr~£  _®2  -  ^  5  ®2  (2 . 5e) 

c 

80 

TjY  =  0  -Bt  <  X  <  -B2  and  B2  <  X  <  Bj  (2.5f) 


For  the  case  of  a  constant  temperature  at  the  lower  surface,  the 
dimensionless  numbers  0  and  Ra^  need  to  be  changed  as 

T  -  Tc 

0  =  f— r  Ra  =  o^3(Th  -  y/01'  (2-5s) 

h  e 

The  corresponding  boundary  condition  (2.5e)  is  replaced  by 

0=1  -B2  <  X  <  B2  (2.5h) 

The  problem  is  specified  mathematically  by  eqns.  (2. 1-2. 5).  The 
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solution  procedure  used  for  solving  these  equations  is  the 
finite- difference  method  SIMPLER,  which  has  been  described  in  detail  by 
Patankar  (1980,  1988).  The  power- law  scheme  is  used  for  the  convection- 
diffusion  terms  and  the  discretization  equations  are  solved  by  using  the 
tridiagonal  matrix  algorithm  (TDMA  or  Thomas  algorithm).  Under- relaxation 
parameters  are  used  to  control  the  advancement  of  the  solutions  and  to 
ensure  convergence.  The  computer  program  was  tested  with  different  grid 
sizes  for  the  same  problem,  and  the  solutions  proved  to  be  essentially 
independent  of  the  grid  sizes. 

2.4  RESULTS  ANT)  DISCUSSION 

2.4.1  Partial  heating  with  constant  heat  flux 

The  numerical  calculations  were  conducted  with  the  configuration  in 

Fig.  2.1,  and  the  results  of  the  temperature  distribution  at  x  =  0  were 

compared  with  the  experimental  data  from  Takenaka  et  al.  (1983).  The 

experiment  was  made  with  a  horizontal  potassium  layer  heated  from  below. 

2 

The  test  vessel  had  a  rectangular  cross  section  of  140  x  96  mm  .  The 

2 

effective  heating  area  in  the  center  of  the  vessel  is  100  x  20  mm  .  The 
vertical  liquid  temperature  distribution  was  measured  by  sliding 
thermocouples  along  the  central  line  of  the  vessel.  Since  the  length  of 
the  heating  element  is  much  larger  than  its  width,  the  heat  transfer  can  be 
modeled  as  two-dimensional  natural  convection  within  the  configuration  as 
shown  in  Fig.  2.1,  with  =  1.70,  Bg  =  0.360,  q  =  10^  W/m^  and  Tc  =  Tg  = 
527  °C.  The  reference  temperature  TQ  is  chosen  as  Tg  and  the  relevant 
properties  are  given  by  Subbotin  et  al.  (1972),  thus  giving  Pr  =  0.004  and 
Ra^  =  3.0  x  10°.  Figure  2.2  shows  the  comparison  of  the  numerical 
temperature  distribution  at  x  =  0  and  the  experimental  data  of  H  =  28  mm 
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Comparison  between  the  numerical  solution 
and  the  experimental  data 


from  Takenaka  et  al.  (1983).  It  can  be  seen  that  the  agreement  is 
generally  good  and  is  excellent  near  the  lower  surface.  The  grid  size  used 
in  the  numerical  calculation  is  33  x  20.  The  temperature  distribution 
consists  of  the  boundary  region  near  the  heating  surface  and  the  bulk 
region.  The  boundary  thickness  is  relatively  larger  compared  with  those  of 
ordinary  fluids.  Figures  2.3  and  2.4  show  the  corresponding  temperature 
contours  and  dimensionless  stream  function  contours,  respectively.  The 
*  isotherms  and  streamlines  indicate  that  a  plumelike  flow  is  generated  above 

the  heated  region,  which  is  the  normal  case  for  natural  convection  in  a 
k  liquid  heated  from  below  with  a  rigid  or  free  upper  surface  (Prasad  and 

Kulacki,  1987;  Boehm  and  Kamyab,  1977). 

Figure  2.5  shows  the  temperature  distributions  for  sodium  with 

different  values  of  the  heat  flux  q.  Since  the  temperature  of  the  upper 

free  surface  is  fixed  at  T£,  higher  values  of  the  heat  flux  results  in  a 

2 

higher  Ttf.  For  example,  the  Tw  of  the  case  with  q  =  30  V/cm  is  almost  50 
o  2 

°C  higher  than  that  of  the  case  with  q  =  5  W/cm  .  Also,  for  the  low  heat 
flux  the  temperature  field  shifts  towards  the  conduction  regime  with  an 
almost  linear  temperature  distribution. 

A  In  the  numerical  calculations,  it  was  found  that  the  convergence  speed 

is  much  slower  with  lower  values  of  the  Prandtl  number,  as  indicated  by 
t  McDonough  and  Catton  (1982)  and  Catton  et  al.  (1974).  Considering  this 

fact,  the  calculations  proceeded  from  higher  to  lower  Prandtl  numbers  with 
the  solutions  of  the  higher  Pr  as  the  starting  values  for  the  lower  Pr. 
Also,  under- relaxation  was  needed  to  ensure  convergent  solutions.  The 
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Figure  2.3.  Calculated  temperature  contours 
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Figure  2.4  Calculated  stream  function  contours 


Temperature  distributions  for  different  heat  flux  q 


under- relaxation  parameter  used  in  the  above  calculations  is  0.1. 


2.4.2  Natural  convection  in  a  horizontal  rigid  cavity  with 
specified  boundary  temperature 

The  problem  of  natural  convection  in  a  rigid  cavity  with  specified 
boundary  temperatures  will  now  be  examined.  The  boundary  conditions  with 
reference  to  Fig.  2.1  in  this  case  are  changed  to  the  following: 
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V  =  U  =  0, 

o 

II 

Y 

II 

1, 

-Bi 

V  =  U  =  0, 

II 
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II 
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II 
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Y  <  1, 

X 

II 

V  =  U  =  0, 
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II 
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Y  <  1, 

II 

X 

<  X  <  B1 

<  X  <  Bx 


The  dimensionless  temperature  and  Raleigh  numbers  are  0  =  (T  -  Tc)/(T^ 
-  T£)  and  Ra  =  g/?H3(Th  -Tc)/ai/,  respectively. 

Raithby  and  Hollands  (1985)  have  proposed  an  empirical  equation  for 
the  natural  convection  heat  transfer  in  horizontal  cavities  nonextensive  in 
the  horizontal  direction,  which  is 

Nu  =  1  +  [1  -  Rac/Ra]*[k1  +  2(Ra1/3/k2)  1_ln  (Ra^3/k2)] 

Ra 

*  [(S380>1/3  '  ']’  (1  -  0.95 [(Ra/Rac)^^  -  1]*})  (2.6) 


where  Rac  is  the  critical  Raleigh  number,  k^  and  k2  are  both  functions  of 
Pr,  and  are  given  as  follows: 
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(2.7) 


kx  =  1. 44/(1  +  0.018/Pr  +  0.00136/Pr2) 

k2  =  75  exp(l .5  Pr'1/2)  (2.8) 

* 

The  square  brackets  with  the  asterisk,  [  ]  ,  indicate  that  only  positive 
values  of  the  argument  are  to  be  taken,  i.e., 

[X]*  =  (|X|  *  X)/2  (2.9)  * 

The  above  equation  has  been  tested  against  experimental  data  for  gases  t 

and  liquids  of  various  Prandtl  numbers  except  liquid  metals,  with  a  maximum 
difference  of  25  percent. 

It  is  of  interest  to  compare  the  present  numerical  solutions  with  eqn. 

(2.6)  and  to  fill  the  gap  in  the  low  Prandtl  number  range.  The 

calculations  were  carried  out  with  Bj  =  0.5  (i.e.,  a  square  box)  and  the 
results  are  presented  in  Fig.  2.6.  It  can  be  seen  that  the  agreement 
between  the  numerical  solutions  and  eqn.  (2.6)  is  generally  good.  The 
maximum  difference  within  the  range  of  comparison  is  less  than  20  percent. 

A  grid  size  of  25  x  25  was  used  in  the  numerical  calculations,  and  the 
solutions  were  obtained  from  higher  Pr  to  lower  Pr  with  the  solutions  of  1 

the  higher  Pr  as  the  initial  guesses  for  the  lower  Pr.  The  number  of 
iterations  for  convergence  ranges  from  700  -  1000.  « 

2.5  CONCLUSIONS 

A  numerical  study  has  been  reported  for  natural  convection  in 
horizontal  liquid  metal  layers  with  localized  heating  from  below.  The 
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2.6),  Raithby  and  Hollands,  (1985 
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Figure  26.  Compaison  of  Eq.  (26)  and  ftie  present  numericd  solutions 


temperature  distribution  of  the  numerical  results  at  x  =  0  agrees  well  with 
the  corresponding  experimental  data  for  a  potassium  layer.  The  isotherms 
and  streamlines  indicate  that  a  plumelike  flow  is  generated  above  the 
heated  region.  The  numerical  results  based  on  a  rigid  square  box  with 
insulated  side  walls,  T^  on  the  bottom,  and  Tc  on  the  top  generally  agree 
with  eqn.  (2.6)  of  Raithby  and  Hollands  (1985),  within  the  Raleigh  number 
range  of  5  x  10^  -  10^,  and  the  Prandtl  number  range  of  0.00125  -  0.01. 


t 
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Section  III 


PERFORMANCE  CHARACTERISTICS  OF  A  THERMAL  ENERGY  STORAGE  MODULE: 
A  TRANSIENT  PCM/FORCED  CONVECTION  CONJUGATE  ANALYSIS 


3.1  SUMMARY 

<  The  performance  of  a  thermal  energy  storage  module  is  simulated 

numerically.  The  change  of  phase  of  the  phase- change  material  (PCM)  and 
4  the  transient  forced  convective  heat  transfer  for  the  transfer  fluid  with 

low  Prandtl  numbers  are  solved  simultaneously  as  a  conjugate  problem.  A 
arametric  study  and  a  system  optimization  are  conducted.  The  numerical 
results  show  that  module  geometry  is  crucial  to  the  design  of  a  space- based 
thermal  energy  storage  system. 


I 


\ 
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3.2  INTRODUCTION 


Recently,  the  study  of  phase- change  thermal  energy  storage  systems  is 
active  due  to  applications  to  space- based  power  plants  and  the  utilization 
of  solar  energy.  Phase- change  materials  (PCM)  have  a  large  latent  heat,  so 
it  is  an  efficient  way  to  absorb  the  heat  energy  during  the  time  period 
when  the  materials  are  subject  to  heat  input  and  to  release  it  to  space 
over  a  long  period  of  time.  A  basic  latent  heat  thermal  storage  geometry 
is  shown  schematically  in  Fig.  3.1.  The  PCM  surrounds  a  pipe,  through  * 

which  the  heat  transfer  fluid  is  passed.  The  problem  by  nature  is  a  time- 
dependent  phase- change  heat  transfer  problem  combined  with  conjugate  forced  * 

convection. 

Heat  transfer  involving  melting  and  solidification  is  a  fertile  area 
for  research  because  of  its  great  importance  in  many  applications.  Since 
problems  of  this  type  are  inherently  nonlinear  due  to  the  existence  of  a 
moving  interface  whose  position  is  not  known  a  priori ,  there  are  relatively 
few  analytical  solutions  to  these  so-called  Stefan  problems.  A  large 
number  of  numerical  techniques  have  been  developed,  but  most  of  the 
numerical  studies  have  focused  on  diffusion- controlled  phase- change 
problems  or  phase- change  problems  including  natural  convection  (Shamsundar 
and  Sparrow,  1975;  Voller  and  Cross,  1981;  Sparrow  and  Ohkuho,  1986;  Voller  t 

and  Prakash,  1987;  Cao  et  al.,  1989a;  Cao  and  Faghri,  1989,  1990). 

* 

The  use  of  a  hollow  cylinder  of  PCM  similar  to  that  in  Fig.  3.1  for  a 
solar  latent  energy  storage  system  was  modeled  by  Solomon  et  al.  (1986).  A 
finite  difference  formulation  with  the  Kirchhoff  temperature  was  used  to 
calculate  the  internal  energy,  temperature,  and  the  position  of  the 
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container  wall 


//// 


PCM 


PCM 


figure  3.1  The  schematic  of  the  latent  heat  storage  system 


phase- change  front.  Stovall  and  Arimilli  (1988)  studied  a  thermal  energy 
storage  system  consisting  of  a  cylindrical  tube  filled  with  a  phase- change 
material  having  a  high  melting  temperature.  The  tube  is  surrounded  by  an 
annular  region  containing  the  liquid  metal  heat  transfer  fluid  for  pulsed 
power  load  applications.  The  emphasis  of  both  of  the  above  studies  is  on 
the  diffusion- controlled  heat  transfer  in  the  PCM.  The  heat  transfer 
between  the  transfer  fluid  and  the  PCM  is  calculated  using  empirical 
correlations  instead  of  solving  the  flow  and  temperature  fields  of  the 
transfer  fluid  numerically  as  a  conjugate  problem.  Also,  the  container 
wall  shown  in  Fig.  3.1  was  ignored.  It  should  be  pointed  out  that  most  of 
the  empirical  correlations  are  based  on  limited  experimental  or  numerical 
data  and  the  fully- developed  conditions;  the  use  of  correlations  may 
increase  the  uncertainty  of  the  numerical  modeling.  Furthermore,  the 
change  of  phase  is  by  nature  a  transient  problem.  The  boundary  temperature 
of  the  transfer  fluid  changes  as  the  phase- change  interface  progresses, 
therefore  the  temperature  field  of  the  transfer  fluid  may  never  reach  the 
fully- developed  state.  This  is  especially  important  for  PCM  storage 
systems  with  liquid  metals  as  the  transfer  fluid.  With  relatively  shorter 
cylinders  and  lower  fluid  velocities,  the  laminar  combined  hydrodynamic  and 
thermal  entry  region  may  dominate  the  flow  along  the  entire  length  of  the 
cylinder. 

In  this  paper,  a  PCM  energy  storage  system  with  the  configuration  in 
Fig.  3.1  is  modeled  numerically.  The  two-dimensional  change  of  phase  for 
the  PCM  and  the  two-dimensional  transient  forced  convection  entrance  region 
for  the  transfer  fluid  with  low  Prandtl  numbers  is  solved  simultaneously. 
Also,  the  influence  of  the  container  wall  is  included  in  the  numerical 
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analysis,  and  an  optimization  of  the  configuration  is  conducted.  To  the 
authors’  knowledge,  this  work  is  the  first  PCM/forced  convection  conjugate 
problem  reported  in  the  literature. 


3.3  Mathematical  Modeling 

The  continuity,  momentum,  and  energy  equations  governing 
two-dimensional  transient  incompressible  laminar  flows  with  no  viscous 
dissipation  in  a  pipe  are  (Ganic  et  al.,  1985) 


1  d  i  dn 

F  It  (rv>  *  S 


(3.1) 


dv  dv  tdv 

at  +  va?  +  ucE 


-Jf  I?  +  vi  tr  I?  (rl?) 


v  5  v-, 

“5  +  7^J 

r  ox 


(3.2) 


du  ir5u  „0u  1  dv  rl  ^  /-fax  ,  Oi 

at  +  va?  +  =  "?f  h  +  vi ar  (ra?)  +  ^ 


(3.3) 


,dT°  ^o  ^o  i  d  ,  dT°,  d2J°, 

pi  cf  <ar  +  var +  uar)  -  kf  t  a?  (rar)  + 


(3-4) 


For  the  PCM,  the  temperature  transforming  model  (Cao  and  Faghri,  1989)  is 
used.  The  energy  equation  is  written  as 


d(Pc°T*)  i  d  n._  w*,  .  d  ,tar\ 
t  =  ?  a?  <kr  ar)  +  a*  (kar) 


at 


(3.5) 


where  T*=  T°  -  T° 
m 
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C°  =  C°(T*)  = 


C_  + 


26T 


S°  =  S°(T*)  = 


cs*T° 

cm  *T°  +  5 


cs  sr  +  h 


(T*  <  -  6T°) 

(solid  phase) 

(-61°  <  T*  <  5T°) 

(mushy  phase) 

(3.6) 

(T*  >  £T°) 

(liquid  phase) 

(T*<  -$T°) 

(solid  phase) 

(-  £T°  <  T*  <  £T°) 

(mushy  phase) 

(3.7) 

(T*  >  0T°) 

(liquid  phase) 

k(T*)  = 


kg  (T*  <  -  £T°) 

ks  +  (k^ks)(T*+5T0)/2W°  (- tfT°  <  T*  <  6T°) 

k ,  (T*  >  tfT°) 


(solid  phase) 
(mushy  phase) 
(liquid  phase) 
(3.8) 


For  the  pipe  wall,  the  energy  equation  is 


n  r  &l°  ,  rl  9  ,  d2T°-i 

pu  cvW  ~  kw  t  TJr  ^  +  ^2  J 


(3.9) 


The  initial  and  boundary  conditions  for  the  case  of  uniform  inlet 
conditions  are  defined  as  follows: 

Initial  conditions:  t  =  0 


Entire  domain:  0  <  x  <  L,  0  <  r  <  r( 
T°  =  T?,  u  =  v  =  0 
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Boundary  conditions:  t  >  0 


Inlet  plane:  x  =  0 


0  <  r  <  r. 
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r  =  r. 


PCM- wall  interface:  0<x<L,  r=r 
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Wall- fluid  interface:  0<x<L,  r  =  r- 
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r  =  r. 
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Outlet  plane:  x  =  L 

•ypO 

°  <  r  <  r»:  g-  -  0 


0  <  r  <  i-.;  =  0 


The  use  of  the  temperature  transforming  model  (Cao  and  Faghri,  1989)  has 
two  advantages.  First,  equations  (3.5  -  3.8)  form  a  set  of  closed- group 
equations,  so  an  explicit,  treatment  of  the  phase- change  interface  is  not 
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needed.  Secondly,  The  time  step  and  grid  size  limitations  are  eliminated, 
which  are  normally  encountered  for  other  fixed- grid  methods.  It  should  be 
pointed  out  that  due  to  the  space  application  the  natural  convection  in  the 
liquid  PCM  has  been  ignored. 


The  following  non-dimensional  variables  are  introduced 


U  D 
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The  dimensionless  continuity,  momentum  and  energy  equations  for  the 
transfer  fluid  are  as  follows: 
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The  dimensionless  energy  equation  for  the  PCM  is 
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1  al  rl  d  /vn57^  ,  d  ,V5I^  dS 
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(3.14) 
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The  dimensionless  energy  equation  for  the  container  wall  is 


51  _  1  flw  rl  d  <?2Tn 

Tt  ~  KS^  LE  M  + 


(3.15) 


The  initial  and  boundary  conditions  are  nondimensionalized  as  follows: 


Initial  conditions:  r  =  0 
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(3.16a) 


Entire  domain:  0  <  X  <  L/D,  0  <  R  <  RQ 
T  =  Tp  U  =  V  =  0 

Boundary  conditions:  r  >  0 


Inlet  plane:  X  =  0 

0  <  R  <  0.5  :  U  =  1,  T  =  1,  V  =  0 
0.5  <  R  <  Rq:  =  0 


Outer  wall:  0  <  X  <  L/D, 

8T 

M 

"o 


PCM- wall  interface:  0  <  X  <  L/D, 


dT 
-  M 


Vail- fluid  interface: 
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0  <  X  <  L/D,  R  =  0.5 


5T 

M 


R  =  0.51 
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r  w 


R  =  0.5' 


Outer  plane:  X  =  L/D 
0  <  R  <  Rq  :  t jjr  -  0 

0  <  R  <  0.5:  =  0 


(3.16b) 


(3.16c) 


(3 . 16d) 


(3.16e) 


(3 . 16f ) 


It  can  be  seen  that  the  temperature  field  can  be  expressed  as 
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* 


♦ 
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T  -  T  (r,  R,  X,  Re^,  Pr^,  o^/a^,  aw/a£j  St,  Cg£,  dT  ,  K^,  kp/ky> 

r„/D,  r0/D,  l/D)  (3.17) 

3.4  Numerical  Procedure 

The  problem  has  been  specified  mathematically  by  eqns.  (3.10  -  3.16). 

The  solution  procedure  used  for  solving  these  equations  is  the 
control- volume  finite- difference  approach  described  by  Patankar  (1980, 
1988).  In  this  methodology,  the  discretization  equations  are  obtained  by 
applying  the  conservation  laws  over  a  finite  size  control  volume 
surrounding  the  grid  node  and  integrating  the  equation  over  the  control 
volume.  The  velocities  and  pressure  are  solved  by  using  the  SIMPLE  scheme 
(Patankar,  1980).  At  the  PCM- wall  and  the  wall- fluid  interfaces,  the 
harmonic  mean  of  the  thermal  conductivity  is  used. 
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The  discretization  equations  are  solved  by  using  the  tridiagonal 
matrix  algorithm  (TDMA  or  Thomas  algorithm).  During  each  time  step, 
iterations  are  needed.  The  converged  results  were  assumed  to  be  reached 
when  the  maximum  relative  change  of  all  variables  between  consecutive 
iterations  was  less  than  0.17..  The  residual  of  the  continuity  equation  was 
also  checked.  The  iteration  was  continued  until  the  sum  of  the  residuals 
was  less  10  .  Different  grid  sizes  for  the  same  problem  have  been  tested 
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and  it  proves  that  both  the  PCM  model  and  the  numerical  scheme  used  for  the 
transfer  fluid  are  essentially  independent  of  grid  sizes  for  the  numerical 
results  in  the  next  section.  The  space  and  time  grid  specifications  are 
given  in  the  next  section  for  the  cases  presented. 

3.5  NUMERICAL  RESULTS  AND  DISCUSSION 

Before  presenting  the  numerical  results  for  the  phase- change  system, 

the  phase- change  model  (eqn.  3.5)  was  checked  against  other  numerical 

results  for  a  two-dimensional  freezing  problem.  Consider  a  liquid 

initially  at  T?  in  a  infinitely  long  prism  with  a  uniform  cross  section 

(Hsiao  and  Chung,  1984).  At  t  >  0,  the  surface  is  kept  at  a  temperature  T° 

w 

<  T°  and  freezing  takes  place  immediately.  Because  of  the  symmetry  of  the 

geometry,  only  a  quarter  of  the  prism  is  considered.  Figure  3.2  shows  the 

o 

interface  position  as  a  function  of  to Ji  along  the  diagonal.  Also 
included  in  this  figure  are  solutions  by  Hsiao  and  Chung  (1984),  and  by  Cao 
et  al.  (1989).  It  can  be  seen  that  the  present  model  agrees  well  with  the 
results  of  these  studies. 

The  numerical  calculations  for  the  thermal  energy  storage  system  were 
then  conducted  with  the  configuration  as  shown  in  Fig.  3.1.  The  system  is 
initially  at  a  temperature  T^  <  0  (less  than  the  melting  temperature  of  the 
PCM).  The  hotter  fluid  enters  into  the  circular  channel  and  heats  the 
system,  which  absorbs  the  energy  from  the  fluid  and  stores  it  as  both 
latent  and  sensible  heat.  The  grid  size  used  for  the  calculation  was  70 
(axial)  x  20  (radial  transfer  fluid)  x  5  (container  wall)  x  15  (PCM)  and  a 
dimensionless  time  step  of  Ar  =  10.  In  order  to  simulate  the  phase- change 
at  a  single  temperature  using  eqn.  (3.14),  the  dimensionless  phase- change 
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temperature  range  £T  is  taken  to  be  0.01. 

Figure  3.3  shows  the  axial  velocity  distribution  in  the  radial 
direction  at  X  =  6  at  different  times.  It  can  be  seen  that  the  velocity 
reaches  the  steady  state  quickly.  After  r  =  45,  the  velocity  profile 
remains  unchanged.  This  does  not  mean  the  velocity  profile  has  reached  the 
fully  developed  condition  along  the  pipe.  Since  the  pipe  is  comparatively 
short  and  the  Prandtl  number  is  small,  the  developing  velocity  region  is  * 

dominant  along  the  whole  pipe  length.  The  temperature  profile,  however,  is 
different.  Figure  3.4  shows  the  radial  temperature  distribution  at  X  =  6 
for  different  time  periods.  The  three  regions  in  the  radial  direction 
(transfer  fluid,  container  wall,  and  the  PCM)  are  also  indicated  in  the 
figure.  The  melting  interfaces  at  different  times  are  the  intersections  of 
T  =  0  and  the  corresponding  temperature  curves.  It  can  be  seen  that  as  the 
melting  interface  progresses,  the  temperature  curve  moves  upward 
accordingly.  Although  the  velocity  field  of  the  transfer  fluid  reaches  the 
steady  state  quickly,  its  temperature  counterpart  cannot  reach  the  steady 
state  before  the  PCM  is  totally  melted.  This  clearly  demonstrates  that  the 
use  of  steady  fully- developed  empirical  heat- transfer  correlations  for  the 
transfer  fluid  may  result  in  a  significant  error  for  the  evaluation  of  the 
system  performance.  * 

Figure  3.5  chows  the  melting  interface  along  the  axial  direction  at  4 

different  times.  At  r  =  1000,  the  melting  interface  has  reached  the  outer 
surface  (r  =  r  )  for  X  <  6,  while  some  PCM  remains  unmelted  for  X  >  6.  The 
reason  is  that  since  the  Prandtl  number  of  the  transfer  fluid  is  very  small 
(the  thermal  conductivity  is  very  large)  a  large  amount  of  heat  is 


60 


H# 


2000 


00  (X)  ^  (N  O  Csl 


O  O  O  O  O 

I 

1 


62 


temperature  distribution  at  X=6  for  different 
time  periods 


o 
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Figure  3.5.  Melting  fronts  along  axial  direction 
for  different  time  periods 


transferred  directly  to  the  PCM  upstream  while  a  relatively  small  amount  of 
heat  is  carried  downstream.  Figure  3.6  shows  the  melting  interfaces  along 
the  axial  direction  for  the  transfer  fluids  with  different  Prandtl  numbers. 
With  a  smaller  Prandtl  number,  heat  transfer  to  the  PCM  is  much  faster  as 
indicated  in  the  figure. 

As  in  normal  forced  convective  heat  transfer,  the  Reynolds  number  has 
a  significant  influence  on  the  system  performance.  This  is  illustrated  in 
Fig.  3.7  with  the  same  Prandtl  number  and  different  Reynolds  numbers. 
Since  the  dimensionless  time  r  contains  the  inlet  velocity  U  ,  the 
comparison  has  been  made  at  a  dimensional  time  t.  The  dimensionless  time 
for  the  case  of  Re^  =  2200  was  taken  as  r  =  1000,  the  corresponding 
dimensionless  times  for  Re^  =  1700  and  Re^  =  1200  were  773  and  545, 
respectively,  with  the  assumption  of  the  same  D  and  for  the  three  cases. 

Since  the  transfer  fluid  is  a  liquid  metal,  the  thermal  conductivity 
has  the  same  order  of  magnitude  as  that  of  the  container  wall.  The 
influence  of  the  thermal  conductivity  of  the  wall  on  the  system  performance 
may  not  be  as  important  as  those  using  gases  as  transfer  fluids,  as 
illustrated  in  Fig.  3.8. 

It  is  very  important  to  evaluate  the  overall  performance  of  the  system 
and  to  optimize  the  geometry  of  the  system  with  given  flow  parameters  and 
thermal  properties.  The  important  system  parameters  used  to  evaluate  a 
thermal  energy  storage  system  are  the  energy  storage  capacity  or  total 
energy  stored,  Q,  (J),  the  energy  storage  density  Q  =  Q./m  (J/kg),  the 
total  latent  energy  stored  (J)  and  the  ratio  of  the  latent  to  the  total 
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Figure  3.8.  Melting  fronts  along  axial  direction  for  different  thermal 

conductivities  of  the  wall  at  r=1000 


1 


energy  stored  Q^/Qt-  The  first  two  parameters,  Qt  and  Qm,  are  most 
important  to  an  energy  storage  system.  In  many  cases,  the  energy  storage 
capacity  is  the  primary  parameter  one  would  be  concerned  with,  while  for  a 
space  application  the  energy  storage  density  is  equally  important  because 
the  weight  of  the  system  is  critical.  With  reference  to  the  geometry  shown 
in  Fig.  3.1,  they  can  be  expressed  as 


fb/D  2  o 
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(3.18) 
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Figures  3.9  and  3.10  present  the  numerical  results  for  the  system 
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Figure  3.9.  Sys 


CM 


o 

mo  ‘*o/'o  *'o  ‘*0 


70 


Figure  3.10.  System  optimization  analysis 
for  different  ro/D  at  t=1000 


optimization  analysis  for  different  L/D  and  rQ/D.  The  dimensional 
parameters  needed  in  eqns.  (3.18  -  3.21)  a-e  H  =  2.9  x  106  J/kg,  p  =  690 
kg/m3,  cp  =  7420  J/(kg  -  K) ,  T°n  -  T°  =  195.4  K  and  D  =  0.1  m.  Both  Qt  and 
follow  the  trend  of  increasing  with  larger  L/D  and  rQ/D,  while  Q^/Qt 
remains  almost  constant.  The  high  proportion  of  latent  heat  storage  in  the 
PCM  is  largely  attributable  to  the  small  difference  between  the  initial 
system  temperature  and  the  PCM  melt  temperature  =  -0.1).  The  trend  of 
the  energy  storage  density  is  opposite  to  that  of  and  Q^.  The  energy 
storage  density  Qm  drops  sharply  with  the  increase  in  L/D  and  rQ/D.  In 
this  situation,  a  trade-off  needs  to  be  reached  when  selecting  the  design 
parameters. 

3.6  CONCLUSIONS 

A  energy  storage  system  with  the  configuration  in  Fig.  3.1  has  been 
studied  numerically.  Numerical  results  show  that  the  fluid  velocity  inside 
the  pipe  reaches  the  steady  state  quickly,  while  the  temperature  field 
continues  to  change  as  the  melting  interface  progresses.  It  is  very 
important  to  treat  the  phase  change  and  the  fluid  flow  as  a  conjugate 
problem  and  solve  them  simultaneously.  The  numerical  results  for  the 
parametric  study  and  geometry  optimization  provide  guidelines  for  the 
design  of  a  space- based  thermal  energy  storage  system. 
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Section  IV 


mathematical  modeling  and  analysis  of  heat  pipe  start- up 

FROM  THE  FROZEN  STATE 


4.1  SUMMARY 

The  start-up  process  of  a  frozen  heat  pipe  is  described  and  a  complete  " 

mathematical  model  for  the  start-up  of  the  frozen  heat  pipe  is  developed 
based  on  the  existing  experimental  data,  which  is  simplified  and  solved  > 

numerically.  The  two-dimensional  transient  model  for  the  wall  and  wick  is 
coupled  with  the  one- dimensional  transient  model  for  the  vapor  flow  when 
vaporization  and  condensation  occur  at  the  interface.  A  parametric  study 
is  performed  to  examine  the  effect  of  the  boundary  specification  at  the 
surface  of  the  outer  wall  on  the  successful  start-up  from  the  frozen  state. 

For  successful  start-up,  the  boundary  specification  at  the  outer  wall 
surface  must  melt  the  working  substance  in  the  condenser  before  dry- out 
takes  place  in  the  evaporator. 


i 
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4,2  INTRODUCTION 


The  heat  pipe  transports  a  large  amount  of  energy  between  a  source  and 
a  sink  by  using  the  latent  heat  of  vaporization  of  the  working  substance 
and  operates  passively  in  the  closed  system  at  various  temperatures 
depending  on  the  working  substance.  The  demand  for  an  effective  thermal 
management  device  for  high  temperature  applications  such  as  cooling  the 
leading  edges  of  reentry  vehicles  and  hypersonic  aircraft,  and  a 
space- based  power  station  stimulates  the  study  of  the  start-up  of  frozen 
liquid- metal  heat  pipes.  Also,  the  start-up  of  frozen  low  temperature  heat 
pipes  is  important  in  applications  such  as  heat  pipe  heat  exchangers, 
cooling  electronic  equipment,  and  melting  the  snow  and  ice  on  roads  and 
bridges. 

Neal  (1967)  and  Shlosinger  (1968)  performed  the  first  experimental 
tests  to  study  the  start-up  performance  of  low  temperature  heat  pipes  with 
the  water  initially  frozen.  Deverall  et  al.  (1970)  also  made  a  series  of 
tests  with  water  and  liquid  metal  heat  pipes.  Successful  start-up  from  the 
frozen  state  was  possible  but  was  highly  dependent  on  the  heat  rejection 
rate  at  the  condenser,  which  should  be  low  enough  to  enable  the  heat  to 
melt  the  working  substance  in  the  condenser  and  allow  liquid  to  return  to 
the  evaporator.  Tolubinsky  et  al.  (1978)  investigated  the  start-up 
characteristics  of  sodium  and  potassium  heat  pipes.  The  characteristics  of 
the  start-up  are  described  based  on  the  experimental  results.  Camarda 
(1977)  investigated  the  performance  of  a  sodium  heat  pipe  cooling  a  leading 
edge.  The  heat  input  was  step-wise  increased  with  time  and  the 
distribution  of  heat  was  non-uniform  to  simulate  aerodynamic  heating.  Heat 
rejection  was  accomplished  by  radiation.  Start-up  and  shut-down  of  a  4  m 
long  lithium  heat  pipe  was  studied  experimentally  by  Merrigan  et  al . 
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(1985,  1986),  in  which  the  heat  input  was  increased  with  time.  A 
temperature  history  of  a  typical  point  on  the  condenser  showed  the  progress 
of  the  effective  heat  pipe  length  with  time.  Ivanovskii  et  al.  (1982) 
presented  the  vapor  temperature  distribution  along  the  length  of  a  sodium 
heat  pipe  during  the  start-up  period.  The  vapor  temperature  was  measured  by 
using  a  movable  microthermocouple  placed  in  the  vapor  channel.  Three  flow 
regimes  in  the  condenser  are  described  based  on  the  vapor  temperature:  free 
molecular  flow,  intermediate,  and  continuum  vapor  flow.  Unfortunately,  the 
existing  experimental  data  for  the  start-up  period  are  in  general  not 
represented  in  the  archival  literature  and  lack  sufficient  information  for 
comparison  with  numerical  simulations.  The  surface  temperature 
distribution  along  the  heat  pipe  length  at  different  times  are  required  in 
addition  to  the  vapor  temperature  to  accurately  analyse  the  start-up 
behavior  of  heat  pipes. 

Colwell  et  al.  (1987)  and  Jang  (1988)  developed  a  simple  mathematical 
model  to  predict  the  start-up  behavior  of  a  sodium  heat  pipe  with  a 
rectangular  cross  section  from  the  frozen  state.  In  the  wall  and  wick 
structure,  energy  transport  is  described  by  the  transient,  two-dimensional 
heat  conduction  equation,  and  the  phase  change  of  the  working  substance  is 
taken  into  account.  In  the  vapor  region,  free  molecular,  choked  and 
continuum  flows  are  considered  and  one- dimensional ,  compressible 
quasi- steady  state  laminar  flow  is  assumed.  The  lengths  of  the  evaporator 
and  condenser  are  not  fixed  and  a  variable  heat  flux  and  radiation  boundary 
condition  are  simultaneously  applied  to  the  outer  surface.  The  numerical 
results  obtained  by  using  the  finite  element  method  are  in  agreement  with 
experimental  results  given  by  Camarda  (1977) . 

The  liquid  metal  heat  pipe  operates  not  only  at  high  temperatures  but 


74 


also  the  initial  temperature  may  be  ambient  temperature.  All  of  the 
experimental  results  mentioned  above  show  that  the  temperature  near  the 
evaporator  increases  rapidly  up  to  a  certain  temperature  and  then  remains 
constant  while  the  temperature  in  the  rest  of  the  heat  pipe  is  unchanged. 
Thus,  a  large  temperature  gradient  exists  between  these  two  regions  and  the 
location  of  this  gradient  moves  toward  the  cooled  end  of  the  heat  pipe  as 
time  increases  until  the  temperature  at  the  end  of  the  condenser  reaches 
the  evaporator  temperature.  In  this  temperature  range,  the  working 
substance  may  be  in  the  solid  state  as  well  as  the  liquid  and  vapor  states. 
In  the  vapor  space,  free  molecular  flow,  continuum  flow,  sonic  and 
supersonic  flow  may  be  encountered  due  to  the  extremely  small  density 
during  the  start-up  of  the  heat  pipe.  These  conditions  may  cause  the 
failure  of  operation  of  the  heat  pipe  and  limit  the  performance  of  the  heat 
pipe.  Understanding  the  start-up  behavior  and  transient  performance  of  the 
high  temperature  heat  pipe  is  therefore  important  and  an  efficient 
mathematical  model  is  needed  to  predict  this  behavior. 

The  first  part  of  this  paper  presents  a  complete  mathematical  model  to 
describe  the  start-up  behavior  of  the  heat  pipe  from  the  frozen  state.  This 
model  is  then  simplified  to  obtain  numerical  results.  The  transient 
two-dimensional  energy  equation  is  solved  numerically  for  the  heat  pipe 
wall  and  wick  structure  saturated  with  the  working  substance  and  the  phase 
change  of  the  working  substance  is  considered.  The  transient 
one- dimensional  compressible  model  is  used  for  the  vapor  flow  dynamics  to 
facilitate  coupling  the  governing  equations  of  the  vapor  with  that  of  the 
wall  and  wick.  To  the  authors’  knowledge,  the  analysis  presented  here  is 
the  only  model  that  includes  the  effect  of  the  transient  vapor  flow  in  the 
analysis  of  the  start-up  of  a  heat  pipe  from  the  frozen  state.  After  the 
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mathematical  models  are  tested  separately,  the  models  are  used  for  a 
parametric  study  of  the  start-up  of  frozen  heat  pipes. 

4.3  DESCRIPTION  OF  EEAT  PIPE  START- PP 

Previous  experimental  observations  suggest  the  following  sequence  of 
events  during  heat  pipe  start-up  from  the  frozen  state.  Initially,  the 
working  substance  is  in  the  solid  state  and  the  vapor  density  is  extremely 
low,  so  that  free  molecular  flow  conditions  prevail  throughout  the  vapor 
space.  The  input  heat  flux  over  the  evaporator  increases  the  temperature 
and  starts  to  melt  the  frozen  substance  in  this  region.  Meanwhile,  the 
heat  transport  from  the  heated  zone  to  the  adjacent  pipe  proceeds  quite 
slowly  via  axial  conduction  through  the  heat  pipe  wall,  working  substance, 
and  wick  structure,  while  the  heat  transfer  in  the  vapor  is  almost 
negligible.  Thus,  a  large  temperature  gradient  exists  between  the 
evaporator  and  condenser. 

When  energy  is  continuously  added  to  the  evaporator,  the  frozen 
working  substance  in  the  evaporator  is  melted,  so  that  evaporation  can  take 
place  at  the  liquid- vapor  interface  and  the  vapor  density  in  this  region  is 
increased.  The  molecular  mean  free  path  in  the  heated  region  then  becomes 
small  compared  to  the  diameter  of  the  vapor  passage  and  the  continuum  flow- 
regime  is  established,  while  in  the  cooled  zone  the  vapor  is  still  in  free 
molecular  flow.  In  the  continuum  flow  region,  the  vapor  flows  into  the 
condenser  section  due  to  the  large  pressure  gradient.  During  this  stage, 
energy  is  mainly  transferred  as  latent  heat  owing  to  vaporization  in  the 
heated  zone,  and  condensation  in  the  cooled  zone  in  the  vapor  space  where 
continuum  flow  is  established.  The  temperature  near  the  evaporator  remains 
constant  and  the  location  of  the  temperature  gradient  moves  toward  the  end 


of  the  condenser  until  continuum  flow  is  establised  in  the  entire  vapor 
space.  Cotter  (1967)  also  described  this  frontal  start-up  mode  when  the 
vapor  density  is  so  low  and  non- condensable  gas  does  not  exist  in  the  vapor 
space. 

In  heat  pipes  with  metallic  working  substances,  the  vapor  densities 
are  very  small  during  the  start-up  even  in  the  continuum  flow  region. 
Thus,  even  for  relatively  low  values  of  heat  input,  sonic  vapor  velocities 
can  be  reached.  Also,  the  vapor  flow  in  the  heat  pipe  is  quite  similar  to 
the  flow  in  a  converging- diverging  nozzle  due  to  the  vapor  addition  in  the 
evaporator  and  the  vapor  removal  in  the  condenser  (Dunn  and  Reay,  1982). 
Thus,  the  heat  transfer  through  the  vapor  space  may  be  limited  by  the 
choked  flow  condition,  and  supersonic  vapor  flow  and  a  shock  front  may 
occur  in  the  continuum  flow  region  in  the  condenser.  The  maximum  rate  of 
heat  transfer  is  limited  by  the  sonic  limit  so  that  a  high  heat  input  in 
the  evaporator  causes  the  various  types  of  start-up  failure. 

This  process  continues  until  the  frozen  working  substance  is 
completely  melted  and  the  continuum  flow  regime  reaches  the  end  of  the  heat 
pipe,  at  which  time  liquid  returning  to  the  evaporator  is  sufficient  for 
normal  transient  operation.  Eventually,  the  heat  pipe  may  reach  a  steady 
state  condition.  The  start-up  process  of  the  liquid  metal  heat  pipe  from  a 
frozen  state  may  be  divided  into  several  distinct  periods  for  convenience 
of  analysis  based  on  the  status  of  the  working  substance  and  the  behavior 
of  the  vapor  flow. 

1.  In  the  first  period,  no  phase  change  takes  place  in  the  entire  region 

bal  the  temperature  near  the  heated  region  increases.  The  vapor  flow 

is  in  the  free  molecular  condition. 

2.  The  working  substance  in  the  evaporator  is  in  the  liquid  state,  but 
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evaporation  does  not  occur  at  the  liquid- vapor  interface. 

3.  The  liquid  and  solid  states  of  the  working  substance  exist 
simultaneously  in  the  wick  structure  and  vaporization  of  the  working 
substance  takes  place  at  the  liquid- vapor  interface.  In  the  vapor 
space,  a  region  of  continuum  flow  is  established  in  the  heated  region 
and  a  continuum  flow  front  moves  toward  the  cooled  end  of  the  heat 
pipe.  The  vapor  flow  may  be  choked  at  the  beginning  of  the  condenser. 

4.  The  working  substance  is  completely  melted  but  free  molecular  flow 
still  exists  in  part  of  the  vapor  space. 

5.  Continuum  flow  exists  over  the  entire  heat  pipe  length  in  the  vapor 
region  but  the  heat  pipe  does  not  reach  the  steady  state  condition. 

6.  The  heat  pipe  then  reaches  the  steady  state  operation. 

For  low  temperature  heat  pipes,  the  experimental  results  of  the 
successful  start-up  from  the  frozen  state  are  very  rare.  Deverall  et  al. 
(1970)  successfully  started  a  water  heat  pipe  from  the  frozen  state  (208 
K) .  The  wall  temperature  distribution  obtained  is  similar  to  that  of  high 
temperature  heat  pipes.  The  vapor  temperature  was  not  obtained,  but  the 
vapor  density  is  relatively  high  even  around  the  melting  temperature.  This 
means  that  the  vapor  velocity  is  very  low  so  that  choked  flow'  and 
supersonic  vapor  velocities  may  not  be  encountered  during  start-up.  Unlike 
high  temperature  heat  pipes,  experimental  results  show  that  the  heat  pipe 
becomes  immediately  active  where  the  ice  is  melted,  but  there  still  is  a 
large  temperature  gradient  in  the  axial  direction.  The  heat' input  at  the 
evaporator  should  be  low  enough  to  return  sufficient  water  into  the 
evaporator  for  the  successful  start-up. 


78 


4.4  MATHEMATICAL  FORMULATION 


Consideration  is  given  to  the  heat  pipe  wall,  the  wick  structuie  with 
the  working  substance  initially  subcooled  at  a  uniform  temperature,  and  the 
vapor  region.  The  external  surface  of  the  heat  pipe  wall  in  the  evaporator 
and  condenser  can  be  exposed  to  heat  flux,  convection,  and  radiation  based 
on  the  application  of  the  heat  pipe.  The  heat  pipe  is  insulated  at  both 
ends  of  the  pipe.  A  schematic  diagram  of  the  physical  model  is  shown  in 
Fig.  4.1. 


4.4.1  Heat  Pipe  Vail 

In  this  region,  the  progress  can  be  modeled  by  the  heat  conduction 
equation  in  a  hollow  cylinder.  The  temperature  and  heat  flux  are 
continuous  at  the  interface  between  the  wall  and  the  wick  structure.  The 
governing  equation  and  boundary  conditions  are  expressed  as  follows: 
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The  initial  condition  is 


Ty(r>X>0)  -  Tq 


(4.2) 


The  applicable  boundary  conditions  at  the  outer  surfaces  of  the  evaporator 
and  condenser  of  the  heat  pipe  are  given  by 
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The  applicable  boundary  conditions  at  the  ends  of  the  heat  pipe  are 


(4.4) 


4.4.2  Vick  Structure  Region 

Initially,  the  working  substance  is  in  the  solid  state  and  the  vapor 
flow  in  the  vapor  space  is  negligible  so  that  the  adiabatic  condition  is 
valid  at  the  liquid- vapor  interface.  When  heat  is  added  to  the  evaporator, 
the  frozen  working  substance  in  the  heated  region  is  melted  so  that  the 
liquid  and  solid  states  of  the  working  substance  exist  in  the  wick.  Fluid 
motion  in  the  liquid  region  may  then  occur  due  to  vaporization  and 
condensation  of  the  working  substance.  The  liquid  flow  in  the  wick  is 
considered  to  be  unsteady  two-dimensional  incompressible  laminar  flow  with 
negligible  body  forces.  The  fluid  and  wick  structure  are  assumed  to  be  in 
local  equilibrium  and  the  velocities  in  the  axial  and  radial  directions  are 
the  local  area- averaged  velocities  over  a  cross  section  of  a  finite  element 
of  the  wick  region  instead  of  the  pore  velocity  or  actual  velocity.  Also, 
the  wick  is  assumed  to  be  isotropic  and  homogeneous.  The  governing 
equations  for  the  wick  region  are  formulated  by  using  the  principles  of  the 
conservation  of  mass,  momentum,  and  energy.  The  continuity  and  momentum 
equations  are  only  effective  when  fluid  motion  exists  in  the  wick.  The 
viscous  dissipation  terms  in  the  energy  equation  are  neglected. 

The  continuity,  momentum  and  energy  equations  are 

i  a 

v  dr  (rV  +  W  =  0  (4*5) 
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where  K  is  the  permeability  of  the  wick  structure,  e  is  the  porosity,  and 
is  the  effective  thermal  conductivity.  The  porosity  is  between  0  and  1 
depending  on  the  porous  material.  When  the  porosity  approaches  unity  (no 
wick  structure  exists),  the  permeability  approaches  infinity.  Therefore, 
eqns.  (4.5- 4.8)  can  be  reduced  to  the  Navier- Stokes  equation  for  unsteady 
two-dimensional  incompressible  laminar  flow.  For  the  steady  state,  eqns. 
(4.5- 4.8)  also  approach  a  special  case  given  by  Hong  et  al.  (1985).  The 
expressions  for  the  porosity  and  the  effective  thermal  conductivity  of 
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the  screen  wicks  are  given  by  Chang  (1987) . 

The  status  of  the  working  substance  during  start-up  is  described 
during  six  distinct  periods  mentioned  in  the  previous  section.  The 
governing  equations  (4. 5-4. 8),  however,  are  not  always  applicable.  For 
example,  during  the  second  period  the  liquid  state  of  the  working  substance 
exists  in  the  wick  but  evaporation  of  the  working  substance  at  the 
interface  is  negligible.  Also,  the  liquid  layer  is  so  thin  that  the  effect 
of  natural  convection  in  the  liquid  region  is  neglected.  Thus,  there  is  no 
fluid  motion  in  the  liquid  region  so  that  only  eqn.  (4.8)  without  the 
second  and  third  terms  is  applicable  and  eqns.  (4.5- 4.7)  are  useful  after 
the  second  period  when  there  is  liquid  motion  in  the  wick. 

In  addition  to  eqns.  (4. 5-4. 8),  coupling  conditions  at  the 

liquid- solid  interface  are  needed: 


Tse  =  he  =  Tm 


(4.10) 
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The  boundary  conditions  at  the  ends  of  the  heat  pipe  are 
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Based  on  the  continuity  of  the  heat  flux  and  temperature  between  the  heat 
pipe  wall  and  the  wick  structure,  the  boundary  conditions  are 
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(4.13) 
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The  boundary  conditions  at  the  liquid- vapor  interface  change  during 
the  start-up  process.  In  the  first  and  second  periods,  free  molecular  flow 
is  prevalent  in  the  vapor  space  so  that  the  boundary  condition  at  the 
liquid- vapor  interface  is 

(4.15) 


In  the  third  period,  the  temperature  of  part  of  the  liquid- vapor 

interface  is  greater  than  the  transition  temperature  from  free  molecular 

* 

flow  to  continuum  flow,  T  ,  and  the  vaporization  and  condensation  of  the 
working  substance  occurs.  New  boundary  conditions  are  then  required  to 
take  into  account  the  phase  change  of  the  working  substance  at  the 
interface.  This  period  will  be  terminated  when  the  frozen  working  substance 
is  completely  melted.  The  new  boundary  conditions  at  the  liquid- vapor 
interface  are  given  as 
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In  the  fourth  and  fifth  periods,  the  working  substance  is  completely 
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melted  so  that  boundary  condition  (4.16)  is  still  valid  at  the  interface 
where  continuum  flow  is  established  and  boundary  condition  (4.17)  is 
expressed  as  follows: 


"t, 

dv 


=  0 


(4.18) 


Vhen  there  is  fluid  motion  in  the  liquid  region,  additional  boundary 
conditions  are  required  for  the  liquid  flow  field.  Boundary  conditions 
(4.13  and  4.14)  at  the  wd.ll- liquid  interface  are  still  effective  and  in 
addition  new  boundary  conditions  resulting  from  the  no- slip  and  impermeable 
conditions  are  given  by 
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In  addition  to  boundary  conditions  (4.16  and  4.17),  the  boundary  conditions 
for  the  velocity  components  at  the  liquid- vapor  interface  where  the 
temperature  is  greater  than  the  transition  temperature  are 


U^(Rv,x,t)  =  Uv(Ry,x,t) 


(4.21) 


5Uv 


(4.22) 
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(4.23) 


Biit  p 

V,(R  j x , t )  = - —  =  V  (R  ,x,t) 

v’  *  ’  pp  Vv  v’  ’  ’  pp 


where  B  is  equal  to  1  for  evaporation,  -1  for  condensation  and  zero  for  the 
adiabatic  section.  The  additional  boundary  conditions  at  both  ends  of  the 
heat  pipe  are  needed  along  with  boundary  condition  (4.12)  as  follows: 


U£(r,0,t)  =  0 

and  U^(r,L,t)  =  0 

(4.24) 

V£(r,0,t)  =  0 

and  V^(r,L,t)  =  0 

(4.25) 

4.4.3  Vapor  flow  dynamics 

Initially,  the  entire  working  substance  is  in  the  solid  state  so  that 
the  vapor  space  may  be  nearly  evacuated.  As  the  temperature  at  the 
interface  increases,  the  vapor  density  also  increases.  Continuum  flow  in 
the  vapor  space  is  considered  to  be  established  when  the  mean  free  path,  A, 
is  substantially  less  than  the  mininum  dimensions  of  the  vapor  flow  passage 
(Holman,  1981),  e.g., 

I„3j<  0.01  (4.26) 

* 

The  transition  temperature,  T  ,  of  the  vapor  corresponding  to  the  given 
dimension  of  the  vapor  space  is  expressed  by  using  the  kinetic  theory  of 
gases  as  follows  (Jang,  1988): 
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(4.27) 


* 

T 


2x10' 4  ^ 


2 


* 

Iterations  are  required  to  obtain  a  value  of  T  due  to  the  temperature 
dependence  of  the  properties.  When  the  vapor  temperature  is  greater  than 
that  calculated  by  eqn.  (4.27),  continuum  flow  is  assumed  to  be  established 
in  the  vapor  space. 

When  continuum  flow  is  established  in  the  vapor  space  during  the 
start-up  from  the  frozen  state,  complex  flow  phenomena  are  encountered  in 
the  continuum  flow  region  due  to  the  extremely  small  density  of  the  vapor. 
The  vapor  pressure  is  low  and  the  temperature  and  pressure  gradients  are 
large  in  the  axial  direction,  so  the  vapor  velocity  may  reach  the  sonic 
velocity,  and  supersonic  vapor  flow  and  a  shock  front  may  occur  in  the 
condenser.  Thus,  the  effects  of  compressibility,  friction  at  the  liquid- 
vapor  interface  and  dissipation  in  the  vapor  should  be  considered  in  the 
mathematical  model.  The  vapor  flow  may  be  considered  to  be  axisymmetric, 
compressible,  unsteady  laminar  flow  and  the  governing  equations  for  this 
flow  are  formulated  with  negligible  body  forces  and  heat  sources  as 
follows.  In  cylindrical  coordinates,  the  continuity,  momentum  and  energy 
equations  are: 


*'v 
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=  0 


(4.28) 
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(4.29) 
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H ^  is  the  second  coefficient  of  viscosity  and  Ey  is  the  total  energy  per 
unit  volume.  The  numerical  and  analytical  solutions  of  the  above  equations 
under  steady  state  conditions  for  annular  and  conventional  heat  pipes  were 
given  by  Faghri  (1986)  and  Faghri  and  Parvani  (1988) .  Transient  results 
are  needed  for  the  start-up  from  the  frozen  condition. 

While  part  of  the  vapor  space  is  in  the  continuum  flow  regime,  free 
molecular  flow  also  exists  in  the  rest  of  the  vapor  space.  Even  though  the 
heat  transfer  through  the  free  molecular  flow  region  may  be  negligible,  the 
boundary  conditions  at  the  border  of  the  two  regions  are  needed  to  solve 
the  governing  equations  for  the  continuum  flow  region.  Since  a  large 
temperature  gradient  exists  in  the  continuum  flow  region  in  the  condenser 
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during  start-up,  most  of  the  vapor  may  be  condensed  at  the  interface.  Thus, 
the  vapor  penetration  in  the  free  molecular  flow  region  may  be  minimal  or 
penetration  may  occur  in  the  immediate  vicinity  of  the  interface  between 
the  continuum  flow  region  and  the  free  molecular  flow  region.  Also,  the 
temperature  in  the  region  of  free  molecular  flow  remains  unchanged  except 
in  the  vicinity  of  the  continuum  flow  region  due  to  the  near  vacuum. 
Therefore,  an  imaginary  plane,  which  is  adiabatic  and  normal  to  the  axial 
direction,  is  assumed  to  divide  the  two  vapor  flow  regions  at  the  point  of 
the  transition  temperature.  The  dividing  plane  moves  towards  the  cooled 
end  of  the  heat  pipe  as  the  location  of  the  transition  temperature  at  the 
liquid- vapor  interface  moves. 

From  the  physical  conditions  at  the  boundary  such  as  the  no- slip 
condition  and  the  adiabatic  condition  at  both  ends  of  the  heat  pipe,  both 
the  axial  and  radial  velocity  components  can  be  assumed  to  be 


uv(r’°’t)  =  U  (r,L,t)  =  V  (r,0,t)  =  V  (r,L,t)  -  0 


(4.32) 


and  the  derivatives  of  the  temperature  are  expressed  as: 

dT  dT 

—jfc  (r>0»t)  =  —ffc  (r,L,t)  =  0  (4.33) 

Also,  at  the  liquid- vapor  interface  of  the  heat  pipe  the  boundary 
conditions  (4.21  -  4.23)  are  effective. 

Since  the  pressure  and  density  at  both  ends  of  the  heat  pipe  are 
unknown,  the  axial  pressure  and  density  gradients  are  used  to  specify  the 
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boundary  conditions  as  follows: 


W 


(r,0,t) 


W 


(r,L,t) 


=  0 


(4.34) 


dp  dp 

(r,0,t)  =  -Tf^-  (r,L,t)  =  0 


(4.35) 


These  conditions  are  realistic  because  the  vapor  velocity  near  the  ends  of 
the  heat  pipe  is  so  low  that  the  variation  of  the  vapor  pressure  is  very 
small . 

Another  boundary  condition  is  induced  from  the  action  of  the  surface 
tension  at  the  liquid- vapor  interface.  The  radius  of  curvature  of  the 
meniscus  is  assumed  to  be  characterized  by  one  radius  of  curvature  at  the 
same  axial  distance.  Then,  the  pressure  difference  at  the  interface  can  be 
given  by  using  the  Laplace  and  Young  equation. 

pv  '  =  ¥  =  r  cos0  (4-36) 

P  c 

Boundary  condition  (4.36)  couples  the  liquid  and  vapor  momentum  equations 
at  the  liquid- vapor  interface.  The  unknown  contact  angle  of  the  liquid  may 
be  evaluated  from  the  conservation  of  the  liquid  mass  in  the  capillary 
structure.  Since  evaporation  and  condensation  of  the  vapor  occurs  at  the 
interface,  the  vapor  temperature  at  the  interface  is  the  saturation 
temperature.  The  Clausius- Clapeyron  relationship  is  used  to  obtain  the 
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saturation  temperature  of  the  vapor  from  the  pressure  as  given  by 


cr 


(4.37) 


The  axisymmetric  condition  along  the  centerline  of  the  heat  pipe 
specifies  the  following  conditions: 


dp  <9U  dE  dp 

'v  V  V 

dr  ~  dr  ~  dr  ~  dr 


(4.38) 


Vv(0,x,t)  =  0 


(4.39) 


4.5  SIMPLIFICATION  OF  THE  MODEL 

The  mathematical  model,  eqns.  (4.1-4.39),  for  the  start-up  behavior  of 
the  liquid  metal  heat  pipe  described  in  the  previous  section  includes  most 
of  the  physical  phenomena  which  may  occur  in  the  heat  pipe.  Therefore,  this 
model  is  very  complex  to  solve  numerically.  The  effect  of  some  physical 
phenomena  may  be  negligible,  so  a  simplified  model  is  derived  to  predict 
the  performance  of  the  heat  pipe  during  the  start-up  period.  For  this 
purpose,  assumptions  are  made  based  on  the  characteristics  of  the  heat  pipe 
and  previous  studies. 

The  density  of  the  liquid  state  of  the  working  substance  is  much 
greater  than  that  of  the  vapor  state,  so  the  velocity  of  the  working 
substance  in  the  wick  structure  is  small.  The  thermal  conductivity  of  the 
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liquid  metal  is  large  and  the  thickness  of  the  wick  region  is  very  thin. 
It  is  then  assumed  that  the  effect  of  the  liquid  flow  in  the  wick  structure 
is  negligible  and  the  wick  structure  is  saturated  by  the  working  substance. 
Thus,  the  heat  transport  through  the  wick  structure  and  working  substance 
is  by  conduction  only  but  the  phase  change  of  the  working  substance  is 
considered.  Under  these  assumptions,  the  same  governing  equation  is  also 
applicable  to  the  heat  pipe  wall  and  the  wick  structure  by  using  the  proper 
properties  for  each  region.  The  governing  equation  is  given  as  follows: 


("Vi 


3T 


d 

3F 


rk. 


Wil 

~5T\ 


d 


k. 

1 


w 


(4.40) 


w  for  wall  region 
se  for  solid  region 
me  for  mushy  region 
le  for  liquid  region 


The  expression 


''  T 


follows: 


("Vi 


<"V» 

'  ("Vfs  *  (1  '  e)  ('Vs 

‘he  {<T  '  V  Mi  -  0  (pcp)s 


l  «  (pcp)f<  +  (1  -  0  0>cp)B 


for  wall  region 

for  solid  region  in 
the  wick 

for  mushy  region  in 
the  wick 

for  liquid  region  in 
the  wick 


(4.41) 
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The  thermal  conductivity,  k^,  can  be  the  thermal  conductivity  of  the  wall 
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material  for  the  wall  region.  When  the  effective  thermal  conductivity  for 
the  wick  region  is  calculated  by  using  the  expression  given  by  Chang 
(1987) ,  the  thermal  conductivity  for  the  working  substance  is  substituted 
by  the  solid  conductivity,  liquid  conductivity,  or  the  average  value  of  the 
solid  and  liquid  conductivity  of  the  working  substance  corresponding  to  the 
solid,  liquid,  or  mushy  region,  respectively. 

The  initial  condition  is 

Tj  (r,x,0)  =  T0  (4.42) 


The  vapor  flow  is  also  simplified  further  from  the  two-dimensional 
model  to  a  one- dimensional  model  since  a  previous  study  (Jang  et  al.,  1989) 
shows  that  the  one- dimensional  transient  compressible  model  described  the 
vapor  flow  dynamics  as  well  as  the  two-dimensional  model  for  simulated  heat 
pipe  vapor  flow.  The  transient  compressible  one- dimensional  continuity, 
momentum  and  energy  equation  are  written  as  follows: 
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(4.44) 


Since  the  governing  equations  are  simplified,  the  corresponding  boundary 
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conditions  are  also  modified  to  match  the  governing  equations.  The 
applicable  boundary  conditions  at  the  outer  surfaces  of  the  evaporator  and 
condenser  of  the  heat  pipe  are  the  same  as  eqn.  (4.3).  The  boundary 
conditions  at  the  ends  of  the  heat  pipe  are 


(4.46) 


The  boundary  condition  at  the  liquid- vapor  interface  where  continuum  flow 
is  not  established  is 


for  T.  <  T* 


(4.47) 


Since  the  one- dimensional  model  is  used  for  the  vapor  flow,  boundary 
condition  (4.16)  at  the  liquid- vapor  interface  is  not  valid.  The  heat 
fluxes,  q(x),  at  the  liquid- vapor  interface  are  calculated  by  using  the 
temperature  distribution  of  the  wick  structure  in  the  radial  direction  and 
then  these  heat  fluxes  are  applied  to  the  one- dimensional  model  of  the 
vapor  flow  as  the  heat  source  and  sink  expressed  by 


q(x)  =  />0VX) 


(4.48) 


The  other  boundary  conditions  for  the  temperature,  velocity,  pressure,  and 
density  at  both  ends  of  the  pipe  for  +he  one- dimensional  vapor  flow  are 
described  by  Jang  et  al.  (1989). 
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4.6  NUMERICAL  PROCEDURES 


The  governing  equation  (4.40)  is  applicable  to  the  heat  pipe  wall  and 
the  wick  structure  and  the  governing  equations  (4.43-4.45)  describe  the 
vapor  flow  dynamics.  For  the  first  and  second  periods,  there  is  no  vapor 
flow  in  the  vapor  space  so  that  only  eqn.  (4.40)  is  effective  with  the 
adiabatic  boundary  condition  at  the  liquid- vapor  interface.  For  the  third 
period,  continuum  flow  is  established  in  part  of  the  vapor  space  so  that 
governing  eqns.  (4.43-4.45)  are  also  solved  and  then  coupled  with  the 
results  for  the  heat  pipe  wall  and  wick  by  using  eqn.  (4.48)  at  the 
liquid- vapor  interface.  Thus,  the  governing  equations  are  separately 
solved  for  each  region.  When  the  coupling  is  implemented  at  the  interface, 
iterations  are  needed.  To  reduce  the  amount  of  computer  time, 
non-iterative  schemes  are  employed  for  each  region. 

The  well-known  alternating  direction  implicit  (ADI)  method  is  used  for 
the  heat  pipe  wall  and  wick,  and  the  phase  change  of  the  working  substance 
during  start-up  is  modeled  by  using  the  equivalent  heat  capacity  method 
(Hsiao,  1985).  This  method  approximates  the  rapid  change  of  the  heat 
capacity  over  the  phase  change  temperature  range,  which  is  an  artificially 
defined  finite  temperature  range,  AT,  instead  of  using  the  Dirac  function. 
Thus,  eqn.  (4.41)  can  be  expressed  as  follows: 
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In  the  numerical  calculation,  this  property  is  evaluated  based  on  the  nodal 
temperatures. 

The  implicit  Beam- Farming  method  is  used  for  the  vapor  flow  dynamics. 
This  scheme  is  a  non-iterative  implicit  method  and  is  similar  to  ADI  for 
multi- dimensional  flow  problems  by  using  factorization.  The  spatial 
derivatives  are  approximated  by  using  the  three- point  second- order  accurate 
central  difference  approximation  for  the  interior  points  and  the  one-sided 
second- order  accurate  difference  approximation  for  the  boundary  nodes.  The 
time  difference  formula  is  second- order  accurate.  This  system  of 
equations  is  solved  using  the  conventional  methods  for  solving  block 
tridiagonal  systems  of  equations.  The  detailed  numerical  method  for  the 
vapor  flow  in  the  heat  pipe  is  described  by  Jang  et  al.  (1989). 

After  continuum  flow  exists  in  the  vapor  space,  eqn.  (4.40)  should  be 
coupled  with  eqns.  (4.43-4.45)  by  using  the  same  boundary  conditions  at  the 
interface.  The  same  amount  of  heat  flux  is  applicable  to  the  governing 
equations  for  the  wick  structure  and  vapor  space  at  the  liquid- vapor 
interface.  Evaporation  and  condensation  occurs  at  this  interface  so  that 
the  temperature  at  the  interface  is  the  same  as  the  saturation  temperature 
of  the  vapor.  Therefore,  the  coupling  of  the  governing  equations  for  the 
vapor  region  to  those  for  the  wall  and  wick  regions  would  be  achieved  by 
using  the  heat  flux  and  the  saturation  temperature  at  the  interface. 
However,  the  heat  flux  at  this  interface  and  the  saturation  temperature  are 
initially  unknown,  so  that  these  boundary  conditions  should  be  assumed  and 
iterations  are  needed  for  each  time  step  until  the  coupling  conditions  are 
satisfied  along  the  interface. 

Initially,  the  heat  flux  at  the  liquid- vapor  interface  is  unknown  when 
continuum  flow  exists.  Therefore,  the  governing  equations  for  the  wall  and 
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wick  regions  are  solved  and  then  the  heat  fluxes  at  each  node  of  the 
interface  are  calculated  based  on  the  temperature  variation  in  the  radial 
direction  in  the  wick  region.  These  heat  fluxes  are  used  as  the  assumed 
boundary  conditions  at  the  interface  for  the  following  time  step.  After 
the  governing  equations  are  solved  for  the  wall,  wick  and  vapor  regions, 
new  heat  fluxes  are  calculated  based  on  the  saturation  temperature  and  the 
temperature  in  the  wick  region.  The  numerical  procedure  used  for  coupling 
is  as  follows: 

1.  It  is  assumed  that  the  liquid- vapor  interface  temperature  is  the 
initial  temperature  for  the  first  time  step. 

2.  Solve  for  the  temperatures  in  the  wall  and  wick  regions. 

3.  Calculate  the  heat  fluxes,  q,  at  each  node  of  the  liquid- vapor 
interface  by  using  the  temperatures  (T.  ^  and  T)  in  the  wick  region. 
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(4.50) 


4.  These  heat  fluxes  are  used  as  the  boundary  conditions  at  the  interface 
to  solve  for  temperatures  in  the  wall  and  wick  regions. 

5.  Use  the  same  heat  fluxes  to  solve  the  vapor  temperature  and  pressure 
for  the  same  period  at  the  wall  and  wick  regions.  Obtain  the 
saturation  temperature,  T  ,  by  using  the  Clausius- Clapeyron 

o 

relationship. 

6.  Calculate  the  new  heat  fluxes,  q',  at  the  interface  by  using  the 
saturation  temperature,  T  ,  in  the  wick. 
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7.  Compare  the  new  heat  fluxes,  q' ,  with  the  old  heat  fluxes,  q,  at  each 
node  of  the  interface. 

8.  If  the  difference  between  the  new  heat  flux  and  old  heat  flux  is 
within  an  acceptable  range,  repeat  steps  4  to  7  for  the  next  time 
step. 

9.  If  the  difference  between  the  new  heat  flux  and  the  old  heat  flux  is 
not  within  the  acceptable  range,  assume  new  guessed  heat  fluxes,  q", 
by  using  the  relaxation  method  and  repeat  steps  4  to  7  until  the 
comparison  of  the  results  is  acceptable. 


q"  =  q  +  a  (q'  -  q)  (4.52) 

10.  Repeat  steps  4  to  9  until  the  temperatures  reach  the  steady  state. 

When  the  coupling  of  the  governing  equations  is  attempted,  some 
physical  characteristics  of  the  two  regions  are  considered.  Since  the 
density  of  the  vapor  is  much  smaller  than  that  of  the  liquid,  the 
volumetric  specific  heat  pc^  (=  8.8  J/m  K)  for  the  vapor  is  much  smaller 
than  that  (=  1040.1  KJ/m^K)  for  the  liquid.  Therefore,  a  difference 
between  the  transient  response  times  of  the  vapor  region,  and  the  wall  and 
wick  regions  exists.  The  time  step  for  the  vapor  space  should  then  be  much 
smaller  than  that  for  the  wall  and  wick  regions.  The  governing  equations 
for  the  wall  and  wick  regions  are  solved  for  one  time  step  by  using  the 
heat  flux  assumed  at  the  interface,  and  then  the  governing  equations  for 
the  vapor  space  are  solved  by  using  the  same  heat  flux  at  the  interface  and 
smaller  time  steps  for  the  same  period  as  the  wall  and  wick  regions.  Also, 
ui1 ike  the  steady  state  case,  the  heat  input  to  the  vapor  space  in  the 
evaporator  is  not  equal  to  the  heat  output  in  the  condenser.  During  one 


time  step  of  the  numerical  procedure,  the  unbalanced  heat  fluxes  are  fixed 
at  the  interface.  The  vapor  temperature  and  pressure  are  very  sensitive  to 
the  heat  flux  due  to  the  small  density  of  the  vapor.  Thus,  the  saturation 
temperature  of  the  vapor  can  be  greater  or  smaller  than  the  wick 
temperature  along  the  axial  length  depending  on  the  difference  of  the  net 
heat  fluxes  at  the  evaporator  and  condenser.  However,  the  saturation 
temperature  is  physically  less  than  the  wick  temperature  in  the  evaporator 
and  is  greater  than  the  wick  temperature  in  the  condenser. 

The  temperatures  at  the  border  between  the  evaporator  and  condenser 

*  change  abruptly  due  to  the  sudden  change  of  the  boundary  conditions  at  the 

external  surface  and  interface.  However,  the  saturation  temperature 
changes  gradually  along  the  axial  distance.  Vhen  the  saturation 

temperature  is  forced  to  become  the  interface  temperature,  a  large 
difference  between  the  saturation  temperature  and  the  wick  temperature  at 
the  interface  exists  at  the  region  separating  the  evaporator  and  condenser. 
Vhen  new  heat  fluxes  at  the  interface  are  calculated  by  using  eqn.  (4.51), 
the  new  heat  fluxes,  q',  for  small  temperature  differences  are  much  greater 
than  the  old  heat  fluxes,  q,  due  to  the  small  distance  (1.2  x  10' ^  m) 

between  nodes  in  the  radial  direction.  This  may  cause  the  numerical 
procedures  to  become  unstable  and  leads  to  the  use  of  small  relaxation 

factors  and  many  iterations. 

»  4.7  RESULTS  AND  DISCUSSION 

The  governing  equations  for  the  wall  and  wick  regions  and  the  vapor 
flow  are  separately  solved,  and  then  are  coupled  at  the  liquid- vapor 

interface  by  using  eqn.  (4.48).  Therefore,  the  numerical  methods  and 
algorithms  can  be  separately  tested  and  compared  to  the  available  data. 
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Since  the  working  substance  changes  phase  from  the  solid  state  to  the 
liquid  state  during  the  start-up  period  from  the  frozen  state,  this  effect 
should  be  incorporated  into  the  numerical  model.  The  solidification  of 
sodium  in  a  square  region  is  chosen  to  verify  the  numerical  model  and 
algorithm  for  governing  eqn.  (4.40).  An  initial  temperature  of  393  K  is 
used,  which  is  greater  than  the  melting  temperature  of  371  K.  A  constant 
boundary  temperature  of  293  K  is  suddenly  imposed  on  the  two  surfaces  and 
the  adiabatic  boundary  condition  is  applied  to  the  other  two  surfaces.  A 
time  step  of  1  second  and  a  uniform  grid  (45  x  45)  is  used  for  the 
numerical  calculation.  The  locations  of  the  interface  between  the  liquid 
and  solid  for  the  time  of  50  seconds  are  calculated  by  interpolating  the 
nodal  temperatures  and  the  results  are  in  agreement  with  available  data 
(Rathjen  and  Jiji,  1971).  The  transient  one- dimensional  model  for  the 
vapor  flow  dynamics  in  the  heat  pipe  has  already  been  verified  by  Jang  et 
al.  (1989).  The  combined  model  is  used  to  predict  the  performance  of  the 
liquid- metal  heat  pipe. 

4.7.1  Simulation  of  the  heat  pipe  start-up  in  the  initial  stages 

Since  the  numerical  scheme  and  algorithm  are  verified  for  the  phase 
change  problem,  the  governing  equations  which  described  the  start-up 
behavior  of  a  cylindrical  heat  pipe  during  the  first  and  second  periods  are 
solved  numerically  for  the  wall  and  wick  regions.  Initially,  the  working 
substance  (sodium)  is  in  the  solid  state  at  the  ambient  temperature.  The 
liquid- vapor  interface  is  assumed  to  be  adiabatic  due  to  the  free  molecular 
flow  in  the  vapor  space  during  these  initial  two  periods.  Thus,  the  heat 
transfer  mainly  takes  place  through  the  wall  and  wick  regions  by  conduction 
only.  The  physical  model  for  the  present  investigation  has  evaporator, 


adiabatic,  and  condenser  section  lengths  of  0.2,  0.1  and  0.2  m, 
respectively.  The  radius  of  the  vapor  space  and  the  inner  and  outer  radii 
of  the  heat  pipe  wall  are  0.00685,  0.008,  and  0.01  m,  respectively.  The 
material  for  the  heat  pipe  wall  and  wick  structure  is  stainless  steel. 

Previous  experimental  results  (Deverall  et  al.,  1970)  show  that 
successful  start-up  of  the  frozen  heat  pipe  greatly  depends  on  the  boundary 
condition  at  the  outer  surface  of  the  evaporator  and  condenser.  The 
working  substance  in  the  evaporator  where  heat  is  added  is  melted  but  that 
in  the  condenser  where  heat  is  rejected  is  still  in  the  solid  state.  Since 
vapor  flow  is  negligible  during  the  initial  stage  of  the  start-up,  a  large 
temperature  gradient  is  developed  between  the  evaporator  and  condenser. 
Also,  the  rate  of  heat  transfer  in  the  axial  direction  is  so  slow  that  the 
working  substance  in  the  condenser  may  not  be  melted  before  available 
liquid  is  vaporized  in  the  evaporator.  Thus,  additional  heat  input  in  the 
evaporator  easily  causes  dry- out  of  the  wick  while  heat  extraction  at  the 
condenser  prevents  raising  the  temperature  of  the  condenser  above  the 
melting  temperature.  Deverall  et  al.  (1970)  made  successful  start-ups  of 
the  frozen  heat  pipe  by  using  a  radiation  boundary  condition  at  the 
condenser.  Since  the  radiation  boundary  automatically  controls  the  heat 
rejection  rate  depending  on  the  temperature  at  the  condenser,  the 
temperature  of  the  condenser  is  allowed  to  rise.  Camarda  (1977)  used  not 
only  a  radiation  boundary  but  also  a  small  amount  of  heat  was  added  on  the 
entire  heat  pipe  surface.  Also,  the  heat  input  was  gradually  increased  as 
time  goes  on.  Unlike  the  transient  behavior  of  the  conventional  heat  pipe, 
start-up  from  the  frozen  state  is  highly  dependent  on  the  state  of  the 
working  substance  in  the  wick.  For  successful  start-up  from  the  frozen 
state,  the  heat  input  and  output  should  melt  the  working  substance  in  the 
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condenser  and  allow  sufficient  liquid  to  return  to  the  evaporator.  Also, 
optimizing  the  start-up  period  is  very  important  to  reduce  possible 
start-up  failure  and  to  perform  the  original  objective. 

In  light  of  the  previous  experimental  results,  the  boundary  conditions 
on  the  outer  surface  of  the  condenser  plays  an  important  role  in  successful 
start-up.  All  of  the  previous  experimental  results  show  the  wall  surface 
temperatures,  so  that  even  for  successful  start-up  the  status  of  the 

working  substance  with  time  is  unknown.  Numerical  simulations  are  ’ 

performed  to  examine  the  effect  of  the  boundary  conditions  and  to  recommend 
the  optimum  boundary  condition.  For  this  purpose,  three  different  boundary  * 

condition  cases  are  chosen  for  the  outer  surface  of  the  heat  pipe  as  shown 
in  Table  4.1.  A  uniform  input  heat  flux  of  50  kV/m  and  radiative  heat 
output  are  used  on  the  evaporator  for  all  three  cases.  An  emissivity  of 
0.9  and  a  radiation  reference  temperature  of  293  K  are  employed.  The 
boundary  conditions  at  the  condenser  and  adiabatic  sections  are  changed 
while  that  at  the  evaporator  remains  the  same  for  all  three  cases.  Sodium 
is  used  as  the  working  substance.  The  melting  temperature  of  sodium  is  371 
K  and  the  transition  temperature  of  680  K  is  obtained  by  using  eqn.  (27) 
for  the  radius  of  the  vapor  space  (0.00685  m)  by  iteration.  For  the 
numerical  calculation,  20  nodes  are  used  at  the  wall  and  wick  regions  in 

the  radial  direction  and  100  nodes  are  used  in  the  axial  direction.  A  time  * 

step  of  1  second  is  used. 

For  case  1,  radiation  is  used  in  the  condenser  to  reject  heat.  Figure  i 

4.2  shows  the  temperature  distributions  at  the  heat  pipe  outer  wall  surface 
and  the  liquid- vapor  interface.  Figure  4.3  shows  the  surface  temperature 
distributions  at  different  times.  As  heat  is  added  in  the  evaporator,  the 
temoerature  in  the  evapoiator  increases  and  at  20  seconds  the  working 
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Table  4.1  Boundary  specifications  at  the  outer 
surface  of  the  heat  pipe  wall 


Case 

No. 

Eva] 

porator 

Adiabatic 

Section 

Condenser 

q(kw/m2) 

Radiation 

Q(kV/n2) 

Radiation 

1 

50.0 

YES 

YES 

0.0 

YES 

2 

50.0 

YES 

YES 

10.0 

YES 

3 

50.0 

YES 

NO 

10.0 

YES 

» 


1 
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Figure  4.2.  Temperature  distribution  at  the  outer  wall  surface  and 
liquid- vappr  interface  of  the  heat  pipe  wall  \yith  tiix^e  for  case  1 
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Figure  4.3.  Temperature  distribution  at  the  outer  wall  surface 
of  the  heat  pipe  wall  with  different  times  for  case  1 


substance  in  the  evaporator  is  in  the  liquid  state.  However  he 
temperature  in  the  condenser  is  not  changed  from  the  initial  temperature 
and  in  the  adiabatic  section  the  temperature  in  the  region  adjacent  to  the 
evaporator  increases  due  to  axial  conduction.  Additional  heat  input  in  the 
evaporator  increases  the  temperature  above  the  transition  temperature  in 
the  evaporator,  but  the  temperature  in  the  condenser  is  still  not  changed 
so  that  a  large  temperature  gradient  exists.  In  the  adiabatic  section, 
part  of  the  working  substance  (sodium)  is  in  the  liquid  state.  When  the 
heat  input  is  continued  at  the  evaporator,  vaporization  occurs  at  the 
interface  in  the  evaporator.  However,  the  working  substance  in  most  of  the 
adiabatic  section  and  condenser  is  in  the  solid  state.  Therefore,  the  heat 
input  in  the  evaporator  should  be  small  to  prevent  dry- out  of  the  wick 
structure  while  the  working  substance  in  the  solid  state  is  melted.  Even 
though  successful  start-up  may  be  possible  for  this  case,  the  start-up 
progresses  very  slowly. 

For  case  2,  10  kV/m  is  added  in  the  condenser  in  addition  to  the 
radiative  boundary  condition  to  assist  in  the  start-up  of  the  frozen  heat 
pipe.  Figure  4.4  shows  the  temperature  distributions  at  the  outer  wall 
surface  and  liquid- vapor  interface.  Figure  4.5  also  shows  the  surface 
temperature  variation  for  different  times.  Since  a  small  amount  of  heat  is 
added  in  the  condenser,  the  temperature  in  the  condenser  is  raised  above 
the  melting  temperature.  However  the  temperature  in  the  adiabatic  section 
is  still  below  the  melting  temperature.  Thus,  liquid  in  the  condenser 
cannot  flow  to  evaporator  until  the  working  substance  in  the  adiabatic 
section  liquefies.  The  temperature  in  the  adiabatic  section  increases 
relatively  faster  than  that  for  case  1  due  to  heat  transfer  at  both  ends  of 
the  adiabatic  section.  The  start-up  period  may  be  shorter  than  that  for 
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Heat  Pipe  Axial  Length,  m 

Figure  4.4.  Temperature  distribution  at  the  outer  wall  surface  and 
liquid- vapor  interface  of  the  heat  pipe  wall  with  time  for  case  2 


Figure  4.5.  Temperature  distribution  at  the  outer  wall  surface 
of  the  heat  pipe  wall  with  different  times  for  case 


case  1. 


Finally,  the  adiabatic  section  is  used  as  part  of  the  condenser  and  10 
kV/m  of  heat  is  input  in  the  condenser  section.  Figure  4.6  shows  the 
temperature  distributions  at  the  heat  pipe  surface  and  the  liquid- vapor 
interface  and  Fig.  4.7  shows  the  surface  temperature  at  different  times. 
Even  though  a  large  temperature  gradient  still  exists  along  the  axial 
direction,  the  working  substance  is  completely  melted  in  the  entire  heat 
pipe  in  80  seconds.  When  vaporization  occurs  in  the  evaporator,  the  working 
substance  can  flow  from  the  condenser  to  the  evaporator  to  prevent  dry-  out 
of  the  wick  structure  in  the  evaporator.  Thus,  a  relatively  large  amount 
heat  can  be  added  at  the  evaporator  without  dry- out  so  that  the  start-up 
period  is  expected  to  be  much  less  than  those  of  cases  1  and  2. 

Based  on  the  numerical  results  for  three  different  cases,  the  boundary 
conditions  for  the  evaporator  and  condenser  should  be  chosen  such  that 
prior  to  vaporization  occurring  in  the  evaporator,  the  working  substance  in 
most  of  the  wick  structure  is  melted  to  be  able  to  provide  liquid. 
Therefore,  a  small  amount  of  heat  input  at  the  condenser  is  recommended  to 
improve  the  start-up  characteristics  and  a  high  heat  extraction  at  the 
condenser  which  prevents  an  increase  in  the  condenser  temperature  should  be 
avoided.  The  presence  of  the  adiabatic  section  also  retards  the  progress 
of  start-up. 

4.7.2  Transient  heat  pipe  operation 

During  the  initial  stages  of  start-up,  the  vapor  space  is  nearly 
evacuated,  so  that  heat  transfer  in  the  vapor  space  is  neglected.  Only  the 
two-dimensional  transient  governing  equation  for  the  wall  and  wick  regions 
are  solved  with  an  adiabatic  boundary  condition  at  the  liquid- vapor 
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Heat  Pipe  Axial  Length,  m 

Temperature  distribution  at  the  outer  wall  surface  and 
or  interface  of  the  heat  pipe  wall  with  time  for  case  3 


interface.  When  continuum  flow  is  developed  in  the  vapor  space,  heat 
transfer  mainly  occurs  through  the  vapor  space  by  using  the  latent  heat  of 
vaporization  and  condensation  at  the  interface.  Thus,  the  one- dimensional 
transient  governing  equation  for  the  vapor  flow  should  be  solved  and 
coupled  with  the  results  for  the  wall  and  wick  regions  at  the  interface. 

However,  the  difference  between  the  densities  of  the  liquid  and  vapor  is  so 
large  that  the  volumetric  specific  heat  of  the  vapor  is  much  less  than  that 
of  the  liquid.  This  difference  leads  to  the  use  of  a  different  time  step 
for  the  transient  governing  equations  for  the  vapor  flow.  Also,  the  vapor 
flow  is  very  sensitive  to  the  heat  input  and  output  at  the  interface. 

During  one  time  step  of  the  numerical  procedure,  the  heat  input  to  the 
vapor  space  is  not  equal  to  the  heat  output  due  to  the  transient  operation 
of  the  heat  pipe.  Thus,  these  factors  make  the  coupling  procedure  difficult 
and  slow. 

To  simulate  the  coupling  of  the  governing  equation  for  the  wall  and 
wick  to  that  for  the  vapor  flow,  the  same  physical  heat  pipe  model  is  used 
except  that  the  adiabatic  section  is  eliminated.  Sodium  is  employed  as  the 
working  substance.  In  order  to  concentrate  on  the  coupling  problem,  it  is 
assumed  that  continuum  flow  is  established  in  the  entire  vapor  space  and 
the  working  substance  is  in  the  liquid  state.  To  yield  these  conditions,  a 
uniform  initial  temperature  of  800  K,  which  is  greater  than  the  transition  ' 

temperature  (680  K) ,  is  used  for  the  wall,  wick,  and  vapor  regions.  The 
external  surfaces  of  the  heat  pipe  wall  at  the  evaporator  and  condenser  are  t 

exposed  to  a  uniform  heat  flux  of  50000  W/m  and  a  convective  boundary 
condition,  respectively.  A  reference  temperature  of  300  K  and  a  heat 
transfer  coefficient  of  100  W/m  K  are  used  for  the  convective  boundary 
condition  at  the  condenser  section.  The  time  step  for  the  vapor  space 
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should  be  much  smaller  than  that  for  the  wall  and  wick  regions,  so  a  time 

step  of  0.1  second  is  employed  for  the  wall  and  wick  regions  and  time  step 
_  ^ 

of  0.1  x  10  c  second  is  used  for  the  vapor  flow.  20  nodes  in  the  radial 
direction  and  160  nodes  in  the  axial  direction  are  used  at  the  wall  and 
wick  regions.  Also,  160  nodes  are  employed  along  the  vapor  space.  A 
relaxation  factor  of  a  =  0.00003  is  used  to  obtain  the  new  guessed  heat 
flux. 

Figure  4.8  shows  the  temperature  distributions  at  the  outer  wall 
surface,  liquid- vapor  interface,  and  the  saturation  temperature  for  a  time 
of  0.3  second.  The  temperature  distribution  at  the  outer  wall  surface  is 
uniform  within  each  section  and  near  the  border  between  the  evaporator  and 
condenser,  the  surface  temperature  abruptly  changes  corresponding  to  the 
boundary  conditions  at  the  surface.  The  interface  temperature  matches  well 
with  the  saturation  temperature  which  is  evaluated  by  using  the 
Clausius- Clapeyron  relationship  with  the  vapor  pressure.  The  saturation 
temperature  decreases  gradually.  Figure  4.9  shows  the  heat  flux 
distribution  at  the  interface.  The  new  heat  flux  calculated  is  converged  to 
the  old  heat  flux.  The  maximum  difference  between  the  two  heat  fluxes  is 
about  107..  Even  though  the  heat  flux  at  the  surface  is  relatively  uniform, 
the  heat  flux  at  the  interface  is  not  uniform.  Also,  the  heat  flux  at  the 
interface  is  much  less  than  that  at  the  surface.  This  implies  that  most  of 
the  energy  is  used  to  raise  the  wall  and  wick  temperature  at  this  moment. 
Figure  4.10  shows  the  vapor  temperature,  pressure,  velocity,  and  density 
distributions.  The  variation  of  the  vapor  temperature,  pressure,  and 
density  is  small.  Also,  a  Mach  number  of  M  =  0.027  is  obtained  at  the  exit 
of  the  evaporator. 
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Temperatures 
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liquid- vapor  interface,  and  saturation  temperature  for  time  of  0.3  second 
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Heat  Pipe  Axial  Length,  m 

Figure  4.9.  Comparison  of  new  heat  flux  with  old  heat  flux  at  the 
interface  for  time  of  0.3  second  during  transient  continuum  flow 


m/dy[  ■  OT  x  ^isuaQ 

St  Vf  L‘Z  ec  6Z  S3 

I _ I _ I _ I _ I _ I 


jaquinN  tpeyj 


SO'T  OO'I  S6'0  06'0  SBO 

2m/N  ‘e_OT  x  9Jnss9Jd 


116 


Heat  Pipe  Axial  Length,  m 

Figure  4.10.  Axial  variations  of  temperature,  pressure,  density  and 
velocity  for  time  of  0.3  second  during  transient  continuum  flow 


4.8  CONCLUSIONS 


The  start-up  process  of  the  frozen  heat  pipe  is  described  based  on  the 
experimental  results.  A  complete  mathematical  model  is  developed  to  predict 
the  start-up  behavior  of  the  heat  pipe  from  the  frozen  condition.  The 
simplified  model  is  used  to  obtain  the  numerical  results.  The  numerical 
results  during  the  first  and  second  periods  show  that  the  heat 
distributions  at  the  outer  surface  of  the  heat  pipe  play  an  important  role 
in  successful  start-up.  The  heat  flux  distributions  for  the  evaporator  and 
condenser  should  be  chosen  to  melt  the  working  substance  in  the  condenser 
prior  to  vaporization  occurring  in  the  evaporator.  A  small  amount  of  heat 
input  at  the  condenser  helps  start-up  and  a  high  heat  rejection  at  the 
condenser  during  the  start-up  should  be  avoided.  The  coupling  of  the 
two-dimensional  transient  model  for  the  wall  and  wick  to  the 
one- dimensional  transient  model  for  the  vapor  flow  is  achieved  when 
continuum  flow  exists  in  the  vapor  space.  During  the  transient  operation, 
the  heat  flux  distribution  at  the  interface  is  quite  different  from  that  at 
the  surface.  Efforts  to  improve  the  present  model  will  continue. 


< 
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Section  V 


A  TEMPERATURE  TRANSFORMING  MODEL  FOR  PHASE- CHANGE 
PROBLEMS  INCLUDING  NATURAL  CONVECTION 


5.1  SUMMARY 

A  temperature  transforming  model  was  proposed  to  convert  the 
enthalpy- based  energy  equation  into  a  nonlinear  equation  with  a  single 
dependent  variable.  The  model  and  the  numerical  scheme  were  first  tested 
with  diffusive  phase- change  problems.  The  model  eliminates  the  time  step 
and  grid  size  limitations  and  is  insensitive  to  the  phase- change 
temperature  ranges,  proving  that  the  model  is  flexible  and  can  handle 
multi- dimensional  phase- change  problems  which  have  large  specific  heat  and 
thermal  conductivity  differences  between  the  liquid  and  solid  phases  in  a 
fixed  grid  system.  The  model  was  further  tested  against  other  phase- change 
problems  including  natural  convection  in  a  cavity.  A  comparison  between 
the  numerical  solutions  and  the  existing  experimental  data  resulted  in  a 
good  agreement. 
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5.2  INTRODUCTION 


Heat  transfer  involving  melting  and  solidification  is  a  fertile  area 
of  research  because  of  its  great  importance  in  applications  including  the 
start  up  of  heat  pipes  from  a  frozen  state,  which  was  one  of  the 
motivations  for  the  present  work.  Since  problems  of  this  type  are 
inherently  nonlinear  due  to  the  existence  of  a  moving  interface,  whose 
position  is  not  known  a  priori,  there  are  relatively  few  analytical 
3  solutions  to  these  so-called  Stefan  problems.  A  large  number  of  numerical 

techniques  have  been  developed  which  can  be  divided  into  two  major  groups. 
t  The  first  group  can  be  called  strong  numerical  solutions  or  classical 

solutions  (Okada,  1984;  Ho  and  Chen,  1986;  Ho  and  Viskanta,  1984;  Sparrow 
and  Ohkubo,  1986).  For  this  group,  deformed  grids  or  transformed 
coordinate  systems  are  required  to  account  for  the  positions  of  the 
phase- change  fronts.  They  are  applicable  to  those  processes  involving  one 
or  two  phases  in  one  space  dimension  which,  with  the  use  of  complicated 
schemes,  are  being  applied  to  two-dimensional  cases  as  well. 

The  second  group  may  be  called  weak  numerical  solutions  or  fixed  grid 
solutions  (Shamsundar  and  Sparrow,  1975;  Raw  and  Schneider,  1985; 

Schneider,  1987;  Voller  and  Cross,  1981;  Voller  et  al.,  1987;  Voller  and 
*  Prakash,  1987;  Prakash  et  al.,  1987;  Cao  et  al.,  1989;  Morgan,  1981; 

Gartling,  1980;  Hsiao  and  Chung,  1984;  Hsiao,  1985).  These  methods  remove 
V  the  need  to  satisfy  conditions  at  the  phase- change  front  and  can  be  easily 

extended  to  multi- dimensional  problems.  In  this  group,  the  two  most 
important  and  widely- used  methods  are  enthalpy  methods  and 
temperature- based  equivalent  heat  capacity  methods.  Both  methods  in  this 
group  have  advantages  and  disadvantages.  Enthalpy  methods  are  flexible  and 
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can  handle  phase- change  problems  occurring  both  at  a  single  temperature  and 
over  a  temperature  range.  The  drawback  of  this  method  is  that  although  the 
predicted  temperature  distributions  and  melting  fronts  are  reasonable,  the 
predicted  time  history  of  the  temperature  at  a  typical  grid  point  may  have 
some  oscillations.  The  temperature- based  fixed  grid  methods  have  no  such 
time  history  problems  and  are  more  convenient  with  conjugate  problems 
involving  an  adjacent  wall,  but  have  to  deal  with  the  severe  nonlinearity 
when  the  phase- change  temperature  range  is  small.  Since  the  objective  of 
this  paper  is  to  improve  the  temperature- based  fixed- grid  methods,  the 
equivalent  heat  capacity  model  is  briefly  described  below.  The  energy 
equation  including  convective  terms  is  given  by  Morgan  (1981)  and  Gartling 
(1980): 


P  Cn  ( 


*  ,  dT°  dT° 


dr 


+  u 


dT 


dT°  0T°  v 

v  W~  w  W  ) 


d 

-fc 


(V 

'  dx  ' 


d 

3y 


dt 


di  (k  w) 


(5.1) 


For  cases  without  convection,  eqn.  (5.1)  reduces  to  the  following  form  as 
used  by  Hsiao  and  Chung  (1984)  and  Hsiao  (1985)  for  phase- change  problems 


*  <?T0  _  d  (,  5T°,  d  (,  5T°  d  ,,  5T°, 
"  W  -  TR  V  W>  *  1?  V  vr>  *  (k  3T1 


(5.2) 


For  the  liquid  phase  and  the  solid  phase,  the  properties  c  and  k  will 
generally  have  different  values.  To  account  for  the  latent  heat  effect  on 
the  liquid- solid  interface,  an  equivalent  heat  capacity  is  introduced 
assuming  that  phase  change  occurs  over  a  temperature  range  as  shown  in  Fig. 
5.1. 
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(T°  <  T°  -  6T°) 

'  m  • 

(T°  .  6To  <  xo  <  To  +  6lo,  (5.3) 

(T°  >  T°  +  6T°) 
v  m  • 

along  with  the  thermal  conductivity  expression 

(T°  <  T°  -  <5T°) 
v  m  ' 

(T°  -61°  <  T°  <  T°  +  6T°)  (5.4) 

(T°  >  T°  +  ^T°) 

Although  the  equivalent  heat  capacity  formulation  has  the  advantage  of 
being  simple  to  program,  many  studies  have  revealed  difficulties  in  the 
selection  of  the  time  step,  mesh  size  and  the  phase- change  temperature 
range,  26T° .  The  predicted  phase- change  interface  location  advances 
frequently  in  an  unrealistic  oscillatory  fashion  and  this  is  accompanied  by 
distortion  of  the  temperature  profile  in  the  region  undergoing  phase  change 
(Hsiao  and  Chung,  1984;  Hsiao,  1985).  Some  attempts  have  been  made  to 
improve  this  method.  For  example,  Hsiao  and  Chung  (1984)  and  Hsiao  (1985) 
used  a  linear  interpolation  of  the  nodal  temperatures  to  account  for  the 
latent  heat  effect.  Although  they  claimed  that  the  algorithm  is 
insensitive  to  the  phase- change  temperature  ranges,  the  time  steps  used  are 
still  small  and  the  numerical  results  are  calculated  with  an  explicit 
scheme  (Hsiao,  1985). 

In  this  paper,  a  new  temperature- based  fixed- grid  formulation  is 
proposed,  and  the  reason  that  the  original  equivalent  heat  capacity  model 
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is  subject  to  such  restrictions  on  the  time  step,  mesh  size  and  the 
phase- change  temperature  range  will  also  be  discussed.  In  addition,  the 
new  methodology  will  be  applied  to  phase- change  problems  including  the 
effect  of  natural  convection  in  a  cavity. 


4 


t 


t 


t 


5.3  NUMERICAL  FORIULATIOW  AND  MODELING 

The  continuity,  momentum,  and  energy  equations  governing 
three-dimensional  transient  laminar  flows  with  no  viscous  dissipation  in 
the  Cartesian  coordinate  system  are  (Kays  and  Crawford,  1980) 


It  +  lx  +  ly  +  dz  ~  0 


(5.5a) 


% t  ^  +  Ec  +  ly  ^uv)  +  !z  ^uw) 


d  i  #ux 


d  ,  chi\  d  ,  5un 

+  +  ^  +  fly  (/*  Jy) 

(5.5b) 


+  M  (puv)  +  ly  + 


d  (  3vv 
+  M  ^  ^ 


dp  d  /  dvs.  d  ,  5vk 

'  +  ph  +  cE  ^  fc)  +  -dy  ^  3y> 

(5.5c) 


+  m  (^uw)  +  ly  ^vw)  +  te  ^w2) 


d  ,  5wx 


d  ,  d\ts.  d  ,  dv\ 
+  p&z  +  W  (**  +  ^ 

(5.5d) 


?(/jE)  5  /  r\  d  i  ^  t  p\  _  d  /i  5T°\  d  /,  5T°\ 

it  +  w  ^uE)  +  ^vE)  +  m  ^wE)  "  m  (k  w)  +  37  (k  vr] 


d  (k 


(5.5e) 


with  the  state  equation 
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(5.6) 


=  c  (T°) 
dT° 

The  momentum  equation  should  be  solved  in  addition  to  the  energy  equation 
even  though  the  focus  in  this  study  is  on  the  energy  equation. 

For  the  phase  change  occurring  over  a  temperature  range  as  shown  in 
Fig.  5.2,  and  in  the  case  of  constant  specific  heats  for  each  phase,  the 
relation  between  enthalpy  and  temperature  can  be  expressed  as 

cs  (T°  '  ”0  +  cs  ^T°  T°  -  T°  <  - 6T°  (solid  phase) 

E(T°)  =  (T°  -  0  +  c.«°  *  \  -n°<-  <  «°( mushy  phase) 

cl  (T°  ■  Tm)  +  CS^T°  +  L  T°  -  T°  >  61°  (liquid  phase) 

(5.7) 

where  T°  is  the  melting  or  freezing  temperature,  and  6T°  =  (T^°  -  Ts°)/2  is 
one- half  of  the  phase- change  temperature  range.  In  the  above  relation,  we 
have  selected  E  =  0  to  correspond  to  phase- change  materials  in  the  solid 
state  at  T°  =  T°  -  £T°.  The  specific  heat  of  the  mushy  phase  has  been 
taken  as  the  average  of  those  of  the  solid  and  liquid  phases;  i.e.,  cm  =  ^ 
(Cs  ♦  c(). 

It  should  be  pointed  out  that  in  the  above  relation  for  the  mushy 
region  a  linear  change  was  assumed.  In  some  real  systems  it  may  take  more 
complicated  forms,  but  this  is  outside  the  scope  of  this  paper. 
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Equation  (5.7)  can  be  rewritten 

with  the 

scaled  temperature 

_  ijiO 

mO 

m 

rc  y  +  cct° 

s  s 

T*  <  -  61° 

(solid  phase) 

E(T*)  = 

(c  +  - — — r)  T*  +  c/T°  +  i 
*  2<5T°  2 

-6T°  <  T*< 

<5T°  (mushy  phase) 

c,  T  +  c  6T°  +  L 
l  s 

T  >  61° 

(liquid  phase) 

(5.8) 


A  temperature  function  is  introduced  as  follows 

E  (T*)  =  C°(T*)  T*  +  S°(T*)  (5.9) 


^  * 

where  C°(T  )  and  S°(T  )  are  determined  from  eqn.  (5.8) 


c  + 

ID 


26  f 


(T*  <  -  6T°) 

(-  «T°  <  T*<  6T°) 
(T*  >  6T°) 


cs6T° 

Cn/T°  + 
cs*T  + 


L 

2 

L 


(T*  <  -  6 T°) 
(-$T°  <  T*<  $T°) 
(T*  >  £T°) 


Upon  substituting  eqn.  (5.9)  into  eqn.  (5.5e) 


(5.10) 


(5.11) 


1  +  H  (puC°T*)  +  ^  (pvC°T*) 

d  /.  3T  a  d  /,  dT  s  „ 

3y  <k  3f)  +  (k  3E~)  +  B 


(5.12) 
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with 


B  =  - 


p(oS°)  »(M°)  .  S(pvS°)  ill vS^l 

dt  dx  +  dy  +  dz 


and  C°  =  C°(T*),  S°  =  S°(T*)  given  by  eqns.  (5.10)  and  (5.11), 
respectively. 


The  energy  equation  has  now  been  transformed  into  a  nonlinear  equation 

* 

with  a  single  dependent  variable  T  .  The  nonlinearity  of  the  phase- change 
problem  is  due  to  the  existence  of  the  moving  phase- change  front.  It  can 
be  shown  that  in  the  liquid  and  solid  phase  regions  away  from  the  moving 
front,  eqn.  (5.12)  reduces  to  a  linear  equation  when  ;he  velocity  field  is 
known . 


It  can  be  seen  that  eqn.  (5.10)  is  identical  to  eqn.  (5.3),  and  that 
when  the  density  p  is  constant  the  equivalent  heat  capacity  model  is  a 
special  case  of  eqn.  (5.12)  with  B  =  0  and  C°  independent  of  the  space 
variables  x,  y,  z  and  time.  This  is  why  many  studies  using  the  equivalent 
capacity  model  found  difficulties  in  selecting  the  time  step  and  mesh  sizes 
and  have  often  encountered  physically  unrealistic  oscillatory  results.  In 

order  for  ^  =  C°  ^  ^ime  step  must 

small  enough  that  the  temperature  change  is  so  small  that  values  of  C°  and 
S°  will  not  change  from  one  phase  to  another  within  the  time  step.  Special 
care  must  also  be  taken  for  the  selection  of  grid  sizes  to  make  C°  and  S° 
independent  of  the  space  variables.  In  fact,  the  above  conditions  are  so 
severe  that  it  is  very  difficult  to  satisfy  them. 
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The  major  difference  between  the  present  model  and  the  equivalent  heat 
capacity  model  is  that  the  present  model  starts  from  the  energy  equation, 
eqn.  (5.5e),  instead  of  eqn.  (5.1).  Since  eqn.  (5.1)  is  only  a  special 
case  of  the  more  general  energy  equation,  eqn.  (5.5e),  the  shortcoming  of 
the  equivalent  heat  capacity  model  is  not  surprising. 


The  conductivity  in  eqn.  (5.12)  is  a  function  of  T  .  With  the 
assumption  of  a  linear  change  in  the  mushy  phase,  k  can  be  written  as 


k(T*) 


(k^  -  ks)  (T*  +  6T°)/2ST0 


(T*  <  -  61°) 

(-61°  <  T*  <  0T°)  (5.13) 

(T*  >  £T°) 


« 


5.4  NUMERICAL  SCHEME  FOR  THE  ENERGY  EQUATION 

Since  the  continuity  and  momentum  equations  are  the  same  as  those  of 
conventional  fluid  flow  problems-  the  emphasis  here  will  be  on  the  energy 
equation.  In  the  following  study,  the  thermophysical  properties  such  as  k 
and  c  are  assumed  to  be  constant  in  each  phase  but  may  differ  among  the 
solid,  mushy  and  liquid  phases,  while  the  density  is  considered  to  be 
constant. 


The  following  non- dimensional ized  variables  are  introduced 


X  - 

T  = 


nO  qO 

CV  o 

-  7 T"  j  3  ”  o  n  5 

c<  - T?) 


ae 

T  -  -5  t  , 

n 


* 
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(5.14) 


„t  '  T?)  ,  cs  ,  ls 

bt - r - ’  Lse  ~  r(  ’  hi  -  iq 


It  should  be  noted  that  the  superscript  o  is  ua.i  for  variables  that  have 
units,  while  unscripted  is  used  for  dimensionless  terms.  The  governing 
equation  is  nondimens ionlized  as 


+  d_  ^UCT^  +  d_  (VCT)  +  d_  (VCT)  =  d_  ^  +  §_  (K  5T^+  d_  ^  0Tj 

Or  d\  dY  dl  dl  dl  dY  dY  dZ  dZ 

+  B*  (5.15) 


with  B*  =  -  f  S  +  —  (US)  +  —  (VS)  +  (VS)  1 
L  0T  ax  a\  ai.  J 


2  (i  +  cs d  + 


2  St  61 


(1  <  -61  ) 

(-61*  <  1  <  61*) 
(T  >  61*) 


(5.16) 


CBltt 


Cs lS1  +  St 


(t  <  -  n  ) 


J  61  (1  +  C  *)  +  mr  (-  <5T  <  T  <  <$T  ) 


(T  >  61  ) 


(5.17) 


«s*  +  ^  -  KsP  (T  +  61 


(T  <  -61  ) 

(-61*  <  T  <  61*)  (5.18) 

(T  >  £T*) 


The  discretization  of  the  above  equation  employs  the  control- volume 
finite- difference  approach  described  by  Patankar  (1980).  In  this 
methodology,  the  discretization  equations  are  obtained  by  applying  the 
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conservation  laws  over  finite  size  control  volumes  surrounding  the  grid 
nodes  and  integrating  the  equation  over  the  control  volume,  which  is  shown 
in  Fig.  5.3.  For  the  grid  point  P,  points  E  and  V  (i.e.,  east  and  west) 
are  its  X-direction  neighbors,  while  points  N  and  S  (i.e.,  north  and  south) 
are  the  Y- direction  neighbors.  The  Z- direction  is  perpendicular  to  the 
paper.  Using  a  fully  implicit  scheme,  adopting  upwind  differencing  in 
space  and  adding  two  more  neighbors  T  and  B  (i.e.,  top  and  bottom)  for  the 
Z- direction,  we  can  write  a  discretization  equation  based  on  eqn.  (5.15) 
for  multi- dimensional  problems 


aPTP  ~  aETE  +  aVTV  +  aNTN  +  aSTS  +  aTTT  +  aBTB  +  b 


where 


(5.19) 


ap  =  Cp  AXAYAZ/Ar  +  Cp  (Fg  +  Fw  +  Fn  +  Fs  +  Ft  +  Fb^  +  De  +  Dw  +  Dn  +  Ds  + 


Dt  +  Db 


aE  "  De  +  CE  Fe 
\  =  Dw  +  CV  Fw 
aN  =  Dn  +  CN  Fn 
aS  =  Ds  +  CS  Fs 
a^,  =  Dt  +  Cj  Ft 

aB  =  Db  +  CB  Fb 


b  =  Cp  Tp  AXAYAZ/Ar  +  (Sp  -  Sp)  AXAYAZ/Ar  -  Sp  (Fg  +  Fw  +  Fn  +  Fs  + 
Ft  +  Fb)  +  SEFe  +  SVFw  +  SNFn  +  SSFs  +  STFt  +  SBFb 

The  flow  rates  and  conductances  are  defined  as 
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K  AYAZ 

F  =  max  [-U  AYAZ,  0] ,  D  =  — - 

6  6  («)e 


K  AYAZ 

F  =  max  [U  AYAZ,  0],  D  =  w— — 

W  W  («)ff 


K  AZAX 

F_  =  max  [-V  AZAX,  0],  D  =  — - 

n  n  n  (SY)n 

K  AZAX 

F  =  max  [V_  AZAX,  0] ,  D  =  -§ - 

s  s  s  (SY)b 


K.  AXAY 

F.  =  max  [-V.  AXAY,  0] ,  D,  =  — - 

t  t  t  (SZ) 


K.  AXAY 

Fb  =  max  [Vb  AXAY,  0] ,  Db  =  JL— 

(M)b 


The  greater  of  a  and  b  is  given  by  max  [a,b] ,  and  Cp,  Tp,  and  Sp 
denote  the  old  values  at  grid  P  for  the  last  time  step.  The  subscripts  E 
and  e,  for  example,  denote  the  values  at  grid  E  and  control- volume  face  e, 
respectively,  for  C,  S,  F,  D,  K,  and  the  velocities. 


Since  eqn.  (5.15)  is  nonlinear,  iterations  are  needed  at  each  time 
step.  This  procedure  is  as  follows: 


f 


1) 


Let  T  represent  the  T  field  as  it  exists  at  the  beginning  of  the  kth 
iteration. 
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2) 


From  these  values,  calculate  tentative  values  of  C,  S  and  K  according 
to  their  relations  with  T,  using  eqns.  (5.16),  (5.17)  and  (5.18), 

respectively. 

3)  Solve  the  nominally  linear  set  of  discretization  equations  to  get  new 

lr-i-1 

values  of  T  . 

4)  Return  to  step  1  and  repeat  the  process  until  further  iterations  cease 

to  produce  any  significant  change  in  the  values  of  T. 

For  the  thermal  conductivity,  K,  it  is  important  to  use  the  harmonic  mean 

described  by  Patankar  (1980)  at  the  faces  of  the  control  volume. 

5.5  NUMERICAL  RESULTS  AND  DISCUSSION 

5.5.1  Phase  change  without  convection 

In  this  case,  the  numerical  scheme  in  the  last  section  is  still 
applicable  with  U  =  V  =  V  =  0.  To  check  the  validity  of  the  present  model, 
the  one- dimensional  and  two-dimensional  numerical  results  are  compared  to 
the  analytical  and  numerical  results  existing  in  the  literature. 

Consider  first  a  one- dimensional  melting  problem  in  a  half- space 
(one- region  problem).  A  solid  at  the  solidification  (or  melting) 
temperature  T°  is  confined  to  a  half- space  x  >  0.  At  time  t  =  0,  the 

temperature  of  the  boundary  surface  at  x  =  0  is  raised  to  T^,  which  is 

higher  than  T°  and  maintained  at  that  temperature  for  t  >  0.  The  exact 
solution  given  by  Ozisik  (1980)  is 

T °(  (x,t)  -  T°  =  (T°  -  I®)  erf  [x/2(o,t)1/2]/erf  (A)  (5.20) 
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for  the  liquid  temperature  distribution,  and 


x(t)  =  2X  (a^t)1/2  (5.21) 

for  the  solid- liquid  interface  position.  The  parameter  X  is  determined  by 
the  relation 


o  i2  c  (Tr  -  T”) 

X2  eA  erf (A)  =  - - - —  (5.22) 

L 

Figure  5.4  shows  the  comparison  between  the  exact  solution  and  the 

present  numerical  solution  with  a  grid  size  AX  =  0.05  for  the  location  of 

the  melting  front  as  a  function  of  the  square  root  of  the  dimensionless 

time  r.  In  order  to  simulate  the  sharp  melting  front,  a  small  temperature 
* 

range  £T  =  0.001  is  used  in  the  calculation  and  the  initial  temperature  is 
set  to  T  =  -0.001  instead  of  T  =  0.  As  can  be  seen  from  Fig.  5.4,  the 
present  solution  is  in  agreement  with  the  exact  solution.  The  error  of  the 
numerical  solution  is  mainly  due  to  the  boundary  condition  of  the  first 
kind  which  gives  a  sudden  temperature  rise  at  r  =  0  and  results  in  a 
nonzero  melting  distance  at  the  onset  of  phase  change.  Also,  with  a  finer 
grid  size,  the  accuracy  will  be  further  increased. 


\ 


t 


The  accuracy  of  the  present  model  will  now  be  tested  with  a 
two-dimensional  freezing  problem.  Consider  a  liquid  initially  at 
temperature  T°  in  a  bar  with  a  uniform  square  cross  section  and  adiabatic 
ends  as  shown  in  Fig.  5.5a.  The  surface  is  suddenly  exposed  to  a  uniform 
wall  temperature  below  the  fusion  temperature  and  freezing  takes  place 
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I  numerical  solutions  of  the  melting 
it  in  a  half-space 


(a)  Bar  of  liquid  with  an  uniform  square  cross  section 


O  T 
1  w 

h~  ^  H 


(b)  One-quarter  of  the  bar  used  for  the  computational 
domain  due  to  the  symmetry  of  the  problem. 

Figure  5.5  Description  of  the  geometry  and  boundary 
conditions  for  the  two-dimensional  freezing  problem 
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immediately.  Because  of  the  symmetry  of  the  geometry,  only  a  quarter  of 
the  bar  is  considered,  as  shown  in  Fig.  5.5b.  The  dimensionless  parameters 
are  k,/ks  =  0.9,  a(/as  =  0.9,  =  (I?  -  T°)/(T°  -  T^)  =  9/7  and  St  = 

(T°  -  T°)/  L  =  2.  Figure  5.6  shows  the  interface  position  as  a  function  of 
dimensionless  time  along  the  diagonal  with  a  grid  size  of  20  x  20  and  a 
time  step  A r  =  0.01.  It  took  about  40  seconds  to  run  the  case  on  a  VAX 
8550.  Also  included  in  Fig.  5.6  are  solutions  obtained  by  Hsiao  and  Chung 

®  (1984),  Keung  (1980),  and  Cao  et  al.  (1989)  using  the  enthalpy¬ 

transforming  model.  It  can  be  seen  that  the  present  two-dimensional 

*  solution  agrees  well  with  the  results  of  those  studies. 

Calculations  were  also  made  with  different  phase- change  temperature 

ranges  and  physical  property  ratios  between  the  solid  and  liquid  phases. 

Figure  5.7  shows  temperature  distributions  along  the  adiabatic  line  of  the 

* 

cross  section  in  Fig.  5.5b  for  different  values  of  AT  .  The  initial 
dimensionless  temperature  of  the  liquid  was  set  to  T^  =  0.5  and  the 
dimensionless  surface  temperature  was  set  to  Ty  =  -0.5.  It  can  be  seen 
that  the  present  scheme  is  insensitive  to  the  phase- change  temperature 
range.  The  temperature  distributions  in  the  solid  region  and  phase- change 
interfaces  (T  =  0)  are  almost  identical  with  slight  differences  of  the 

*  temperature  distribution  in  the  liquid  region.  It  can  also  be  seen  that 

*  * 
the  curve  with  AT  =  0.01  is  almost  the  same  as  that  with  AT  =  0.001. 

*  This  means  that  a  moderately  small  phase- change  temperature  range  is  enough 
to  simulate  the  phase- change  problem  occurring  at  a  single  temperature. 

Time  step  and  grid  size  limits  were  not  encountered  n  the  above 
calculations.  The  dimensionless  time  step  A r  can  be  on  the  order  of  0.1  - 
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=  9/7 


Figure  5.6  Interface  position  along  the  diagonal 
with  prescribed  boundary  temperature 


Temperature  profiles  along  an  adiabatic  line  for 
different  values  of  6T*  and  grid  sizes 


1.0  and  the  real  time  step  At  can  be  as  large  as  several  days.  The  model 

is  also  independent  of  grid  sizes.  The  curve  with  a  grid  size  of  40  x  40 

* 

is  almost  identical  to  that  of  20  x  20  with  the  same  AT  .  The  CPU  time  for 

* 

running  the  case  with  AT  =  0.1  in  Fig.  5.7  up  to  r  =  2.0  was  less  than  17 

seconds  on  the  VAX  8550  with  a  time  step  of  A r  =  0.2  and  a  grid  size  of 

* 

2u  *  20.  The  CPU  time  for  running  the  case  with  AT  =  0.001  in  Fig.  5.7  up 
to  r  =  2.0  was  about  5  minutes  with  a  time  step  of  Ar  =  0.2  and  a  grid  size 
of  40  x  40  on  the  same  computer.  The  scheme  can  deal  with  phase- change 
problems  with  large  specific  heat  and  thermal  conductivity  differences 

between  the  solid  and  the  liquid  phases,  which  is  often  the  case  in  many 

applications. 

It  should  br  pointed  out  that  due  to  the  severe  nonlinearity  of 

phase- changes  problems,  underrelaxation  is  needed  during  the  iterations  for 

* 

C,  S  and  K  in  eqn.  (5.15)  when  AT  is  very  small.  The  normal 

* 

underrelaxation  parameters  are  0.1  -  0.5  for  small  values  of  AT  .  For  a 

* 

moderate  value  of  AT  ,  underrelaxation  is  not  needed. 


5.5.2  Phase  change  with  natural  convection 

To  further  demonstrate  the  present  model,  the  phase  change  with 
natural  convection  in  a  two-dimensional  cavity  as  shown  in  Fig.  5.8  was 
studied.  The  dimensionless  governing  equations  with  parameters  introduced 
in  eqn.  (5.14)  and  the  Boussinesq  approximation  are  as  follows 


AU  +  AV  = 
AX  AY 


(5.23) 
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(5.24) 


*5  +  «fl?l  +  ism  =  .  »  +  L  (Pr  «S)  +  L  (Pr  «V) 

dr  dl  dY  dl  dl  dl  dY  6Y 

«v  +  im  +  djfl  __ .  m  +  Ra  p  T  t  L  {Pr  «,  +  L  (Pr  £)  (5.25) 

dr  5X  dY  dY  *  dX  dX  dY  dY 

^11  +  fL  (UCT)  +  —  (VCT)  =  —  (K  — )  +  —  (K  — ) 

dr  dl  dY  dl  dl  dY  dY 

-  01  +  L  (US)  +  —  (VS)1  (5.26) 

Ldr  ax  ay  -• 

where  Ra  =  g/?H^  (T^  -  T°)/i/^  a^,  Pr^  =  z^/a^,  Pr  =  v/a and  C,  S,  and  K 
are  given  by  eqns.  (5.16),  (5.17)  and  (5.18),  respectively.  Equations  (5.23 
-  5.26)  were  solved  numerically  using  the  iterative  SIMPLER  algorithm 
(Patankar,  1980)  with  the  numerical  scheme  in  the  previous  section  of  this 
paper  used  for  the  energy  equation.  The  SIMPLER  algorithm  is  based  on  a 
fully- implicit  discretization  scheme  for  the  unsteady  terms.  The 
discretized  form  of  the  momentum  equation  is  very  similar  to  eqn.  (5.19). 
An  important  difference  is  that  the  grids  used  are  "staggered"  over  the 

temperature  grid.  A  consequence  of  the  staggered  grid  approach  is  that 

care  has  to  be  taken  in  numerically  implementing  momentum  sources  and 
physical  properties  which  depend  on  temperature. 


At  the  solid- liquid  interface  or  in  the  mushy  region,  the  velocities  U 
and  V  go  to  zero.  If  the  momentum  equations  above  are  to  be  applicable  in 
phase- change  problems,  this  behavior  must  be  taken  into  account.  There  are 
a  number  of  alternative  methods  by  which  the  behavior  of  the  velocities  in 
the  mushy  region  can  be  modeled  (Voller  et  al.,  1987).  A  commonly  used 
procedure  is  to  prescribe  a  fluid  viscosity  which  is  equal  to  +be  liquid 


♦ 


* 
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viscosity  in  the  liquid  region  and  which  increases  gradually  over  the  mushy 
zone  to  a  large  value  in  the  solid.  This  representation  can  be  used  quite 
conveniently  and  was  adopted  in  the  present  study  because  finite- difference 
procedures  can  properly  handle  large  discontinuities  in  the  diffusion 
coefficients.  Therefore,  the  viscosity  is  expressed  as 

vi  (T  >  ST) 

u£  +  (T  -  ST*)  (i/£  -  N) /2ST*  (-ST*  <  T  <  ST*)  (5.27) 

N  (T  <  -  ST*) 

where  N  takes  a  sufficiently  large  value  such  as  10^.  Notice  that  the 

term  Ra  Pr^T  in  eqn.  (5.25)  is  independent  of  v.  Prior  to  the  calculation 

of  the  phase- change  problems,  the  computer  program  written  with  the  SIMPLER 

algorithm  was  tested  with  pure  natural  convection  in  the  cavity,  and  the 

results  with  a  grid  size  of  25  x  25  were  compared  with  those  given  by 

Raithby  and  Hollands  (1985)  with  an  agreement  of  5  percent  within  the  whole 

3  5 

given  Raleigh  number  range  (10  -  10  ) . 


t 


T 


Calculations  were  conducted  with  a  one- region  melting  problem 

occurring  at  a  single  temperature  (Okada,  1984).  The  initial  and  boundary 

conditions  for  this  problem  are  T^  =  0,  T^  =  1,  T£  =  0.  In  order  to 

simulate  the  phase  change  having  a  sharp  interface  with  the  present  model, 
* 

£T  was  taken  as  0.01  and  T^  and  Tc  were  set  to  -0.01  in  the  calculations. 
Figure  5.9  shows  the  calculated  results  with  a  20  x  20  grid  size  compared 
with  the  experimental  results  given  by  Okada  (1984).  Also  included  in  the 
figure  are  the  numerical  results  with  the  coordinate  transforming  method 
and  the  quasi- steady  assumption.  It  can  be  seen  that  the  present  results 
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Comparison  of  the  locations  of  the  melting  front 


are  in  agreement  with  the  past  experimental  and  numerical  results.  The  CPU 
time  for  running  the  case  is  about  3  hours  on  the  VAX  8550,  which  is  much 
higher  than  those  of  pure  diffusive  phase  change  problems. 

The  one- region  restriction  is  relaxed  and  c  and  k  were  permitted  to 
take  different  values  between  the  liquid  and  solid  phases.  The  calculation 
was  conducted  again  with  T^  =  0.5,  T^  =  0.5,  Tc  =  -0.5,  Ra  =  10^,  Pr^  =  30, 

T  St  =  0.5,  Cg^  =  0.5,  Kg^  =  1.5  ,  £T  =0.01,  and  a  grid  size  of  30  x  30. 

The  numerical  results  for  this  two- region  freezing  problem  at  r  =  2.0  are 
i  presented  in  Fig.  5.10  and  Fig.  5.11  by  a  dimensionless  temperature  contour 

and  a  stream  function  contour,  respectively. 

The  advantage  of  the  present  scheme  is  that  the  calculations  can  be 
made  in  a  fixed  grid  system  without  involving  complicated  coordinate 
transforming  procedures  and  many  assumptions.  The  present  scheme  can  also 
be  applied  to  three-dimensional  problems. 

5.6  CONCLUSIONS 

The  temperature  transforming  model  and  the  numerical  scheme  proposed 
in  this  paper  proves  to  be  flexible,  easy  to  implement,  and  able  to  handle 

*  various  complicated  phase- change  problems  including  those  with  large 

specific  heat  and  thermal  conductivity  differences  between  the  solid  and 

*  liquid  phases.  For  diffusion- controlled  phase- change  problems,  the  model 

eliminates  the  time  step  and  grid  size  limitations  which  the  equivalent 
heat  capacity  model  usually  encounters  and  is  insensitive  to  phase- change 
temperature  ranges.  The  model  can  also  properly  deal  with 

convective- diffusive  phase- change  problems  in  a  fixed  grid  system.  The 
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present  solutions  have  been  compared  with  previous  analytical  and  numerical 
results  as  well  as  the  experimental  data  with  a  good  agreement.  The 
advantage  of  this  model  is  that  explicit  attention  does  not  have  to  be 
given  to  the  nature  of  the  phase- change  front  and  the  model  can  be  applied 
to  three-dimensional  problems  without  involving  cumbersome  mathematical 
schemes . 
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Section  VI 


A  TRANSIENT  TVO- DIMENSIONAL  COMPRESSIBLE  ANALYSIS  FOR  HIGH  TEMPERATURE 
HEAT  PIPES  WITH  A  PULSED  HEAT  INPUT 


6.1  SUMMARY 

The  transient  behavior  of  two-dimensional  heat  pipes  is  studied 
numerically.  The  energy  balance  relation  for  the  liquid- vapor  interface  is 
derived  and  a  simplification  concerning  the  porous  wick  is  justified  by  a 
simple  analysis.  The  numerical  results  for  compressible  vapor  flow  with 
high  Mach  numbers  are  compared  with  the  experimental  data  in  the 
literature.  The  transient  responses  of  heat  pipes  to  a  pulsed  heat  input 
are  then  elaborated.  The  numerical  results  show  that  for  transient  heat 
pipe  analysis,  it  is  very  important  to  include  the  porous  wick  and  the  wall 
in  the  numerical  calculations,  and  to  treat  the  entire  heat  pipe  as  a 
single  system  rather  than  analyze  the  vapor  flow  alone. 
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6.2  INTRODUCTION 


The  study  of  heat  pipe  dynamics  is  important  for  start-up,  shut-down 
and  operational  transients  such  as  pulsed  heat  loads.  Most  of  the  previous 
works  on  the  transient  behavior  of  heat  pipes  treat  the  vapor  flow  with 
one- dimensional  models  which  incorporate  friction  coefficients  found  from 
two-dimensional  numerical  results  or  experiments  (Ivanovskii  et  al.,  1982; 
Jang  et  al.,  1989).  Heat  pipes  are,  in  general,  intended  for  isothermal 
operation  at  a  higher  vapor  temperature,  but  for  the  transient  operation, 
the  two-dimensional  character  of  heat  pipes  is  evident.  Also,  the  friction 
coefficient  correlations  used  in  the  one- dimension; 1  vapor  models  are  based 
on  limited  experimental  or  numerical  data  and  their  validity  is  not  fully 
confirmed.  Therefore,  the  use  of  empirical  correlations  will  surely 
increase  the  uncertainty  of  the  models. 

Bowman  and  Hitchcock  (1988)  and  Bowman  (1987)  studied  the  transient 
two-dimensional  compressible  flow  of  air  in  a  simulated  heat  pipe.  The 
experimental  data  for  the  simulated  steady  state  heat  pipe  vapor  flow  were 
obtained  and  compared  with  the  numerical  results  with  a  fairly  good 
agreement.  Both  the  experimental  data  and  the  numerical  scheme  are 
valuable  to  understand  the  flow  vapor  structure  and  pressure  drops  at  high 
Mach  numbers.  However,  the  vapor  flow  was  not  actual  vapor  flow  in  a  heat 
pipe  but  air  flow  in  a  porous  pipe.  Both  the  heat  pipe  wick  and  wall  were 
also  not  considered  in  the  transient  numerical  simulation.  Because  the 
time  period  of  the  vapor  transient  is  very  short  when  compared  to  that  of 
the  heat  pipe  wick  and  wall,  the  numerical  results  can  hardly  simulate  the 
actual  transient  processes  in  a  heat  pipe. 

Issacci  et  al.  (1988)  studied  the  transient  behavior  of  the  vapor  flow 
in  a  rectangular  heat  pipe  with  the  SIMPLE  numerical  scheme  developed  by 
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Patankar  (1980) .  Flow  reversal  was  detected  in  the  condenser  and  adiabatic 
regions  with  water  as  the  working  fluid.  However,  only  some  axial  and 
radial  velocity  profiles  were  presented  in  the  paper,  and  no  information 
concerning  vapor  compressibility  and  the  transient  behavior  of  a  heat  pipe 
was  given.  Also,  the  heat  pipe  wall  and  wick  were  not  included  in  the 
numerical  analysis. 

In  this  paper,  a  liquid- vapor  interface  relation  will  be  derived,  and 
the  transient  and  steady  two-dimensional  compressible  vapor  flow  including 
the  effects  of  the  wick  and  wall  will  be  studied  numerically. 

6.3  IATHEIATICAL  FORMULATION  AND  MODELING 

The  heat  pipe  configuration  and  coordinate  system  studied  in  this 
paper  is  presented  in  Fig.  6.1.  The  wick  is  saturated  with  the  liquid 
phase  of  the  working  fluid  and  the  remaining  volume  of  the  tube  contains 
the  vapor  phase.  Heat  applied  at  the  evaporator  causes  the  liquid  to 
vaporize  into  the  vapor  space.  The  vapor  flows  to  the  condenser  and 
releases  the  latent  heat  as  it  condenses.  The  released  heat  is  rejected 
into  the  environment  by  convection  or  radiation  from  the  outer  condenser 
surface. 


6.3.1  Governing  equations 

The  governing  equations  for  transient  compressible  laminar  vapor  flow 
with  constant  viscosity  are  as  follows  (Ganic  et  al.,  1985): 

It  +  F  I?  ^rv)  +  !i  =  0  (6A) 
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Figure  6.1  The  heat  pipe  configuration  and  coordinate  system 
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The  equation  of  state  is  given  by 


p  =  pRT 


(6.5) 


More  effort  is  needed  to  describe  the  boundary  conditions  for  the 
*  vapor  flow.  In  fact,  the  most  challenging  problems  arise  from  the 

treatment  of  the  liquid- vapor  interface  and  the  liquid  flow  in  the  porous 
,  wick.  The  governing  equations  for  the  heat  conduction  in  the  wick  and  the 

pipe  wall  are  also  presented  in  the  following  section. 


6.3.2  Analysis  of  the  vapor- liquid  interface 

Figure  6.2  shows  a  control  volume  surrounding  the  vapor- liquid 
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Liquid 


I 


Vapor 


control  volume 


Figure  6.2  Vapor-liquid  interface  control  volume 


interface.  The  thickness  of  the  control  volume  should  be  considered  to  be 
infinitesimal.  In  the  figure,  is  the  heat  flux  from  the  liquid  side. 
The  convective  heat  flux  due  to  the  vapor  velocity  at  the  interface  is 


‘conv 


.  c  T  =  p-v-c  T 
1  pv  s  ^1  1  pv 


pv  s 


The  heat  flux  due  to  heat  conduction  into  the  vapor  is 
qcond  v  dn 

1  where  n  is  the  normal  direction  of  the  interface.  The  latent  heat  due  to 

phase  change  from  the  liquid  phase  to  the  vapor  phase  or  vice  versa  is 

«i  =  "i  hfg  =  f i  vi  hfg 


where  Tg  is  the  saturation  temperature  at  the  interface,  v^  is  the  radial 
vapor  velocity  at  the  interface,  p ^  is  the  interfacial  vapor  density,  and 
nu  is  the  corresponding  mass  flux. 


t 


v^>0  (condensation) 
v^<0  (evaporation) 

0  (adiabatic) 


Making  an  energy  balance  at  the  interface,  we  get 


qin  qconv  +  qcond  +  qL 
91 

vi  =  (qin  +  kv  W)  /  (hfg  +  cpv  Ts)  pi 


16.6) 
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o 

Ve  have  neglected  the  kinetic  energy  v^/2  in  the  above  energy  balance 
because  the  velocity  at  the  interface  is  comparatively  small. 

The  problem  remains  how  to  express  q^n,  which  is  composed  of  two 
parts: 

^in  ”  qcond,£  +  ^conv,£  "  "^eff  cln-  +  p t  vt  cp^  ^s 

The  first  term  on  the  right-hand  side  in  the  above  equation  is  the 
heat  conduction  from  the  liquid  side  into  the  interface  control  volume,  the 
second  term  is  the  convective  heat  flux  due  to  the  liquid  velocity  v^, 
where  k^^  is  the  effective  conductivity  in  the  porous  liquid- wick,  which 
can  be  calculated  with  appropriate  relations  given  by  Dunn  and  Reay  (1982). 

According  to  the  conservation  of  mass,  we  have 


n  vi  =  pi  vi, 


vt = 


n Vi 


Since  p ^  «  p it  follows  that  v^  «  v^ 

It  should  be  pointed  out  that  because  v^  is  very  small,  its  influence 
on  the  vapor  velocities  at  the  interface  boundary  is  negligible.  However, 
the  term  p^  v^  may  not  be  small  because  p^  is  large. 

Let  us  now  evaluate  the  effect  of  the  liquid  mass  flow  on  the 
temperature  distribution  in  the  liquid  wick.  Since  the  porous  wick  is  vei y 
thin  (on  the  order  of  1  mm)  and  the  thermal  conductivity  of  liquid  metals 
is  high,  the  temperature  distribution  in  the  wick  could  be  assumed  lirear 
as  shown  in  Fig.  6.3,  where  T^w  is  the  temperature  at  the  wick- wall 
interface,  6^  is  the  thickness  of  the  wick,  and  Tg  is  the  saturation 
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Figure  6.3  Temperature  distribution  in  the  wick 
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temperature  at  the  vapor- liquid  interface. 
For  steady  state  operation,  we  have 


^cond,£  ^in  ^v£^scp£  ^eff  5^ 


T.  -T 
1W  s 


I  =!inWtft,  , 
lw  Vf  1  s 


Now,  suppose  qin  is  only  attributed  to  heat  conduction  in  the  wick, 


i.e. , 


91 1 

^in  ~  ^eff  chi 


Then,  the  temperature  distribution  in  the  wick  must  be  somewhat 
different.  Let  T-  be  the  wick- wall  interface  temperature  in  this 

1W 

situation: 


Tiw  =  £ 


in 


eff 


VTS 


Therefore,  the  maximum  temperature  difference  is 


T.  =  "JlMi  s 

1K  keff  ' 


For  a  typical  heat  pipe,  6^  =  0.5  mm,  vi  =  5  m/s,  Tg  =  818  K,  = 

1262  J/(kg-K),  pi  =  0.461  x  10' 2  kg/m3,  kgff  =  35  V/(m-K). 
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“CftP  s 
keff 


p-v-c  tT 

s.  =  PJ^  el 5 

1  *eff 


^  =  0.34  K 


It  can  be  seen  that  the  maximum  temperature  difference  between  these 
two  cases  is  negligible.  This  is  why  many  studies  found  that  the 
conduction  model  is  appropriate  to  describe  the  heat  transfer  in  the  wick 
for  liquid  metal  heat  pipes  (Chen  and  Faghri,  1989).  It  is  also  worth 
pointing  out  that  the  above  conclusion  is  only  applicable  to  the  working 
fluids  that  have  a  high  thermal  conductivity.  For  working  fluids  with  low 
thermal  conductivities  such  as  water,  the  conduction  model  may  need  some 

* 

modifications. 

Based  on  the  above  analysis,  we  will  neglect  the  liquid  velocity  in 

the  wick  and  apply  the  conduction  model  to  the  heat  transfer  in  the  wick. 

This  does  not  mean  that  the  liquid  flow  in  the  porous  wick  is  not 

important.  In  fact,  the  liquid  flow  in  the  porous  wick  is  very  important 

to  determine  the  capillary  limit  of  heat  pipes.  Since  the  major  concern  in 

this  paper  is  the  heat  pipe  dynamics,  the  porous  wick  is  assumed  to  be  so 

designed  that  it  has  enough  capillary  force  to  drive  condensate  to  the 

evaporator  as  long  as  the  working  fluid  in  the  wick  is  melted.  Therefore, 

for  the  porous  wick  and  wall,  the  energy  equation  (6.4)  is  still  applicable 
DP 

with  v=w=0,j^r=f=0  and  the  corresponding  properties  Cp,  p,  k  of  the 
wick  and  wall,  i.e. , 

^Veff  M  =  keff  ^F  If  +  (for  the  wick)  (6-7) 

pv  Cpw  If  =  kw  t?  If  (f!f)  +  ^  (for  the  wal1)  ( 6 -8) 
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where  (/>cp)eff  =  uPic^l  +  (1"40  (/,cp)m>  u  porosity  of  the  wick,  and 

(/>Cp)m  is  the  heat  capacity  of  the  wick  structure  material. 

6.3.3  Boundary  conditions 

At  both  ends  of  the  heat  pipe,  the  no- slip  condition  for  velocity  and 
the  adiabatic  condition  for  temperature  are  applied. 

z  =  0  and  L: 


v  =  w  =  0 


(6.9) 


At  the  centerline,  the  symmetry  condition  implies 


« 


r  =  0: 


(6.10) 

For  the  boundary  conditions  at  the  vapor- liquid  interface  (r  =  Ry) , 
the  evaporator,  adiabatic  and  condenser  sections  are  described  separately. 
For  the  evaporator  and  adiabatic  sections,  the  temperature  at  the 
vapor- liquid  interface  is  assumed  to  be  the  saturation  temperature 
corresponding  to  the  interface  pressure: 


3w 

3F 


=  0,  v  =  0 


5T 

3r 


=  0 
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(6.11) 


where  pQ  and  TQ  are  the  reference  saturation  temperature  and  pressure. 
The  interface  velocities  following  the  analysis  are 


r  =  Ry  (evaporator  and  adiabatic  sections) 

91,  91 

vi  “  (~keff  W~  +  kv  W)  /  (hfg  +  cpvTs)  pi  (6'12) 

w  =  0 


It  is  worth  noting  that  because  of  the  axial  conduction  in  the  wall 
and  wick,  the  evaporator  and  adiabatic  sections  are  not  necessarily  the 
same  as  those  prescribed  at  the  outer  wall  surface.  The  actual  location  of 
the  adiabatic  section  at  the  interface  is  where  the  values  of  v^  calculated 
with  the  equation  above  approaches  zero. 

For  the  condenser,  the  interface  velocity  is  usually  described  as  a 
mass  suction  velocity,  but  this  is  not  exactly  true.  The  interface 
velocity  is  not  only  determined  by  the  exterior  factors  but  also  by  the 
*  interior  flow  field  of  the  vapor.  Because  of  the  closed  fluid  circulation 

system  of  the  heat  pipe,  the  mass  condensed  at  the  condenser  interface  must 
*  be  equal  to  that  evaporated  at  the  evaporator  interface  for  the  steady 

state  operation.  This  is  somewhat  similar  to  the  outflow  boundary 
condition  of  a  normal  duct  flow.  Therefore,  the  upwind  scheme  is  used  and 
the  diffusion  is  neglected  at  the  liquid- vapor  interface  of  the  condenser 
section.  For  the  velocity  component  w,  the  no- slip  boundary  condition  is 
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still  applicable,  w  =  0,  at  r  =  Ry. 

At  the  condenser  interface,  vapor  condenses  and  releases  its  latent 
heat  energy.  The  latent  heat  energy,  as  well  as  the  convective  heat 
energy,  is  absorbed  at  the  interface.  In  order  to  simulate  this  process,  a 
heat  source 


«s  =  <hfg  *  cpv  Tiv>  »i  vi  (6-13> 

is  applied  at  the  interface  grids,  where  T^v  is  the  vapor  temperature  at 
the  grids  next  to  the  interface  on  the  vapor  side,  and  v.  can  be  obtained 
by  a  mass  balance  over  the  grid. 

Since  the  vapor  flow  is  solved  by  using  the  SIMPLE  algorithm  (this 
will  be  described  in  the  next  section),  the  normal  velocities  at  the 
boundaries  are  given  and  the  staggered  grid  arrangement  is  used,  so  no 
information  concerning  p  at  the  boundaries  will  be  needed  (Patankar,  1980). 
The  boundary  condition  at  the  liquid- wall  interface  is 


r 


k  51  -  k  91 
w  dr  eff  dr 


(6.14) 


For  the  outer  pipe  wall  surface: 


Qrn 

evaporator  kw  ^ 


r  =  R  ^e 
o 


adiabatic 


ffl 

di 


r  =  R 


=  0 


\ 

(6.15) 


(6.16) 
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r  h  (T  -  T  )  (convection)  (6.17a) 

C  V  a 

r  =  Rq  <r  eT^  (radiation)  (6.17b) 

where  h  is  the  convective  heat  transfer  coefficient,  and  T  and  T  are  the 

outer  wall  surface  temperature  and  the  environment  temperature, 

respectively,  e  is  the  emissivity  and  a  is  the  Stefan- Boltzmann  constant. 

It  should  be  pointed  out  that  for  the  heat  pipe  dynamics,  T  should  not  be 

w 

t  fixed,  it  will  be  allowed  to  change  during  the  transient  process. 

,  6.4  NUMERICAL  PROCEDURE 

The  governing  equations  along  with  the  boundary  conditions  were  solved 
by  employing  the  control- volume  finite- difference  approach  described  by 
Patankar  (1980,  1988).  The  vapor  flow  inside  the  heat  pipe  was  solved 
using  the  SIMPLE  algorithm.  There  are  two  alternatives  to  deal  with  the 
compressibility  of  the  vapor  flow  (Patankar,  1980;  Karki,  1986).  One 
method  is  to  use  the  density- correction  formulation,  in  which  the  vapor 
density  is  expressed  as 

p  -  p  +  p'  =  p  +  Kp'  (6.18) 

<  where  K  =  dp/d p,  p  =  p/RT  and  the  subsequent  coefficient  of  the  pressure 

correction  equation  consists  of  two  terms,  i.e.,  a  diffusive  term  and  a 
»  convective  term  due  to  density  correction.  The  second  method  is  to  choose 

the  pressure  as  a  dependent  variable  and  directly  apply  the  state  equation 
p  =  p/RT  to  obtain  the  vapor  density  while  iterating.  Based  on  our 
numerical  results,  it  was  found  that  the  latter  is  more  convenient  and  the 
speed  of  convergence  is  faster  when  the  Mach  number  is  high.  With  moderate 


condenser  ,  Sfl 

~  w 
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Mach  numbers  (Ma  <  1),  both  methods  give  comparable  results.  Therefore, 
the  second  method  was  selected  in  our  numerical  calculations.  The  sequence 
of  numerical  steps  is  as  follows: 


1.  Initialize  the  pressure,  temperature  and  velocity  fields.  The  density 
values  are  obtained  from  the  current  pressure  and  temperature  fields 
through  the  state  equation. 

*  * 

2.  Solve  the  momentum  equations  to  obtain  v  and  w.  1 

3.  Solve  the  p'  equation  and  update  the  current  pressure  field. 

4.  Calculate  v,  w  from  their  starred  values  using  the  velocity- correction  « 

formulas . 

5.  Solve  the  temperature  equation. 

6.  Steps  (1)  through  (6)  are  repeated  until  convergence  is  reached. 

The  discretization  equations  for  w,  v,  T  and  p'  have  the  general  form 

a?(J>p  =  ag(j)g  +  ay<J)y  +  a^<j)j^  +  ag(J)g  Jr  b  (6.19) 

which  is  described  by  Patankar  (1980) . 

Since  the  values  of  the  density  will  normally  be  available  only  at  the 
main  grid  points,  the  interface  densities  need  to  be  calculated  by  a  * 

closure  relation.  In  this  case,  an  upwind- biased  density  is  employed  as 

suggested  by  Karki  (1986).  The  density  at  interface  e,  for  instance,  is  * 

calculated  using  the  following  expression 

Pe  =  Pp  max  [ve/|we|,  0]  +  max  [-  wg/ 1 wg |  ,0]  (6.20) 
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A  combination  of  the  direct  method  (TDMA)  and  the  Gauss- Seidel  method 
was  employed  to  solve  the  discretization  equations.  The  converged  solution 
was  assumed  to  be  reached  when  the  relative  change  of  all  variables  was 
less  than  10’^. 

6.5  NUMERICAL  RESULTS  AND  DISCUSSION 

6.5.1  Compressible  vapor  flow  in  a  simulated  heat  pipe 

Attention  will  first  focus  on  the  simulated  heat  pipe  vapor  flow,  and 
the  numerical  results  will  be  compared  with  the  experimental  data  of  Bowman 
(1987)  to  verify  the  mathematical  model  and  algorithm.  The  experimental 
data  was  obtained  by  simulating  the  vapor  flow  of  a  cylindrical  heat  pipe 
with  a  porous  pipe  which  has  an  inside  diameter  of  1.65  m  and  a  length  of 
0.61  m.  The  evaporator  and  condenser  sections  have  equal  lengths  and  were 
simulated  by  the  injection  and  suction  of  air  without  phase  change  at  the 
interface.  Since  the  experiment  was  focused  on  the  vapor  flow,  the  porous 
wick  and  heat  pipe  wall  were  not  considered. 

Figure  6.4  shows  the  numerical  results  of  the  steady  axial  pressure 
profiles  along  the  centerline  compared  with  the  corresponding  experimental 
data  by  Bowman  (1987).  It  can  be  seen  that  the  agreement  between  the 
experimental  data  and  the  numerical  solutions  is  generally  good.  Figure 
6.5  shows  the  corresponding  axial  Mach  number  variations  along  the 
centerline  of  the  simulated  heat  pipe.  For  the  two  subsonic  cases,  the 
axial  velocity  increased  until  it  decelerated  in  the  condenser  due  to  mass 
removal.  For  the  supersonic  case,  mass  removal  caused  a  further 
acceleration  of  the  air  in  the  condenser  until  a  shock  was  encountered. 
The  pressure  profiles  are  similar  to  those  in  a  converging- diverging 
nozzle.  Kemme  (1969)  and  Bystrov  and  Popov  (1978)  reported  observing 
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Figure  6.4  The  axial  pressure  profile  along  the  centerline  of 
the  simulated  heat  pipe  for  steady  state 


ixial  Mach  number  profile  along  the  centerline  of 
simulated  heat  pipe  for  steady  state 


a  similar  phenomenon  by  measuring  the  axial  temperature  variation  along  a 

sodium  heat  pipe.  It  should  be  pointed  out  that  for  Ma  >  1,  considerable 

axial  temperature  and  pressure  changes  will  exist  and  the  heat  pipe 

operation  will  be  far  from  normal  working  conditions.  This  is  indeed  the 

reason  that  most  heat  pipe  designs  use  Ma=l  as  the  sonic  limit  of  heat 

pipes.  Since  in  the  experiment  the  uniform  source  and  sink  pressures  for 

the  blowing  and  suction  regions  were  specified,  the  radial  mass  flow  rate 

along  the  pipe  at  the  interface  is  not  uniform  as  indicated  by  Bowman  ' 

(1987),  so  a  linearized  radial  mass  flow  rate  has  been  specified  at  the 

interface  to  match  the  experimental  boundary  conditions.  « 

The  above  calculation  has  been  made  with  a  grid  size  of  20  x  60. 

Different  grid  sizes  have  been  used  to  calculate  the  same  problem  and  it  is 
found  that  the  program  is  essentially  independent  of  grid  sizes  for  the 
above  specifications.  The  steady  state  was  assumed  to  be  reached  when  the 
relative  change  of  all  variables  was  less  than  0.1  percent.  The  CPU  time 
for  running  a  case  from  an  arbitrarily  guessed  initial  condition  to  the 
steady  state  is  about  5-26  minutes  on  a  VAX  6320. 

6.5.2  Transient  response  of  the  heat  pipe  to  a  pulsed  heat  input 

In  this  section,  a  heat  pipe  with  configuration  shown  in  Fig.  6.1  is 

studied.  The  heat  pipe  is  treated  as  a  single  system,  and  the  vapor  flow,  * 

porous  wick  and  the  heat  pipe  wall  are  solved  as  a  conjugated  problem.  The 
heat  pipe  has  a  stainless  steel  wall  and  sodium  is  the  working  fluid.  The  • 

dimensions  and  properties  are  L  =  0.105  m,  L  =  0.0525  m,  L  =  0.5425  m, 

Ry  =  0.007  m,  6^  =  6y  =  0.001  m,  ky  =  0.0352  V/m-K,  k  ff  =  45  V/m-K,  kw  = 

21.7  V/m-K,  cpv  =  2580  J/kg-K,  cp^  =  1262  J/kg-K,  and  hf  =  4.227  x  106 
J/kg.  The  outer  wall  boundary  condition  at  the  condenser  was  first  chosen 
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as  a  convective  boundary  condition  (eqn.  6.17a),  with  the  heat  transfer 

o 

coefficient  h  =  39.0  M/m  -K  and  an  environment  temperature  T  =  300  K. 

C  d 

The  calculation  was  conducted  with  a  heat  input  Q  =  623  V  at  the  outer  wall 
of  the  evaporator.  The  steady  state  vapor  temperature  at  the  centerline  is 
shown  in  Fig.  6.6  and  is  labeled  as  t  =  0.0.  At  this  time,  a  sudden 
increase  in  heat  input  from  Q  =  623  V  to  Q  =  770  V  was  imposed  at  the 
evaporator.  The  transient  respAuse  of  the  vapor  temperature  is  presented 
*  in  the  same  figure.  Since  the  heat  transfer  coefficient  at  the  outer  wall 

of  the  condenser  is  constant,  with  a  larger  heat  input  Q,  the  outer  wall 
»  temperature  at  the  condenser  needs  to  be  raised  to  reject  more  heat  into 

the  environment.  Since  the  outer  wall  temperature  is  directly  related  to 
the  vapor  temperature  in  the  heat  pipe,  the  vapor  temperature  will  increase 
accordingly.  After  about  10  minutes,  the  heat  pipe  reached  another  steady 
state  as  shown  in  the  figure.  It  is  noticeable  that  at  a  higher  working 
temperature,  the  vapor  temperature  is  more  uniform  along  the  heat  pipe  than 
that  at  a  lower  working  temperature.  The  time  step  used  in  the  above 
calculation  was  At  =  0.5.  It  took  about  57  minutes  of  CPU  time  on  the  CRAY 
401R  with  a  grid  size  of  20  x  40. 

In  space  applications,  radiation  boundary  conditions  are  more  often 
encountered.  The  calculation  was  conducted  again  for  the  same  heat  pipe 
<  with  a  radiative  boundary  condition  (eqn.  17b,  e  =  0.85)  at  the  outer  wall 

surface.  The  transient  vapor  temperature  is  presented  in  Fig.  6.7.  Like 
»  the  case  with  the  convective  boundary  condition,  the  curve  labeled  with  t  = 

0.0  is  the  steady  state  vapor- temperature  with  Q  =  623  V.  At  t  =  0,  the 
b^at  input  is  suddenly  increased  to  Q  =  770  V  and  remains  constant 
afterwards.  In  this  case,  the  transient  time  is  shorter  and  the  heat  pipe 
reaches  another  steady  state  in  about  6  minutes.  Figure  6.8  shows  the 
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temperature  along  the  centerline  for  different  time 
with  the  convective  boundary  condition 


temperature  along  the  centerline  for  different  time 
with  the  radiative  boundary  condition 
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corresponding  outer  wall  temperature  along  the  heat  pipe.  Although  the 
vapor  temperature  tends  to  be  more  uniform  at  a  higher  working  temperature, 
its  outer  wall  counterpart  is  still  stepwise.  This  trend  is  reasonable  and 
agrees  with  experimental  observations.  As  will  be  illustrated  in  Fig. 
6.13,  since  heat  needs  to  be  conducted  to  the  vapor- liquid  interface 
through  the  pipe  wall  and  the  wick,  the  outer  wall  surface  temperature  at 
the  evaporator  is  higher  than  the  vapor  temperature,  while  at  the  condenser 
heat  needs  to  be  conducted  from  the  liquid- vapor  interface  to  the  outer 
surface  to  dissipate  into  space,  so  the  vapor  temperature  at  the 
interface  is  higher  than  that  of  the  outer  wall  surface. 

As  the  heat  pipe  approaches  the  steady  state,  the  heat  rejected  at  the 
outer  wall  surface  of  the  condenser  must  approach  the  input  heat  Q  at  the 
evaporator.  This  energy  balance  has  been  checked  with  the  numerical 

rLc 

results.  At  the  steady  state,  the  difference  between  Q  and  Q  =  2ir 

c  o  J0 

e  a  T^  dz  is  less  than  17.. 

w 

The  corresponding  transient  response  of  the  density,  pressure  and 
axial  velocity  for  the  heat  pipe  with  the  radiation  boundary  conditions  are 
presented  in  Figs.  6.9-6.11,  respectively.  It  can  be  seen  that  both 
density  and  pressure  increase  with  the  pulsed  heat  load,  while  the  axial 
velocity  behaves  differently.  Shortly  after  the  pulsed  load  is  imposed, 
the  axial  velocity  increased  slightly  and  then  decreased  gradually  to 
»  another  steady  state  condition.  Since  the  vapor  is  compressible,  m  = 

Vj^A,  the  rate  at  which  the  density  increased,  at  some  time,  exceeded  that 
of  mass  evaporation  at  the  evaporator;  therefore,  the  vapor  velocity  began 
to  decrease  as  shown  in  Fig.  6.11. 

Figures  6.12  and  6.13  presented  the  radial  distributions  of  axial 
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pressure  profile  along  the  centerline  of  the  heat  pipe 
for  different  time  periods 


?nsity  along  the  centerline  of  the  heat 
different  time  periods 
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Figure  6.1Z  Axid  velocity  distributions  in  the  radid  direction  at  different  axid 

locations  and  times 


velocity  and  temperature,  respectively.  Different  regions,  namely,  the 
vapor  space,  the  porous  wick,  and  the  heat  pipe  wall  are  also  indicated  in 
the  Fig.  6.13.  As  discussed  concerning  the  outer  wall  temperature 
distributions,  the  temperature  slope  is  positive  in  the  evaporator, 
negative  in  the  condenser,  and  zero  in  the  adiabatic  section.  Also,  at  a 
higher  working  temperature,  the  vapor  temperature  tends  to  be  more  uniform 
in  the  radial  direction,  especially  at  the  evaporator.  It  should  be 
pointed  out  that  the  heat  transfer  mechanisms  are  different  in  different 
regions.  In  the  regions  of  the  heat  pipe  wall  and  wick,  heat  is 
transferred  mainly  by  diffusion,  while  in  the  vapor  region,  heat  is 
transferred  mainly  by  convection. 

From  the  numerical  results  presented  above,  it  is  clear  that  both  the 
transient  response  and  the  steady  state  of  a  heat  pipe  heavily  depend  on 
the  boundary  conditions  at  the  outer  wall  surface  of  the  condenser.  This 
is  illustrated  by  comparing  the  numerical  results  in  Figs.  6.6  and  6.7, 
where  the  only  difference  is  the  boundary  conditions  on  the  condenser.  The 
temperature  distributions  indicate  that  the  extent  of  the  departure  from 
uniform  temperature  along  the  heat  pipe  is  determined  not  only  by  the  heat 
load,  but  also  by  the  working  temperature  level  which  is  dependent  on  the 
boundary  condition  on  the  condenser. 

Many  of  the  previous  numerical  analyses  on  heat  pipe  dynamics  were 
focused  on  the  vapor  flow  alone,  and  it  was  found  that  the  vapor  transient 
periods  were  very  short.  For  example,  the  vapor  transient  time  without 
considering  the  wick  and  wall  as  reported  by  Bowman  and  Hitchcock  (1988) 
and  by  Issacci  et  al.,  (1988)  is  only  about  0.001  ~  0.01  seconds.  It  seems 
plausible  that  a  steady- state  vapor  model  could  be  used  in  conjunction  with 
a  transient  wall  model  to  study  the  heat  pipe  dynamics,  but  this  is  not  the 


179 


case.  The  numerical  results  in  this  paper  show  that  the  vapor  transient 
period  have  the  same  magnitude  of  5  ~  10  minutes  as  that  of  the  heat  pipe 
wall.  This  is  because  the  vapor  and  the  wall  are  coupled  through 
liquid- vapor  interface  by  evaporation  and  condensation.  Therefore,  it  is 
very  important  to  solve  them  simultaneously. 

The  advantage  of  the  present  numerical  modeling  is  that  it  is  mainly 
based  on  the  energy  balance,  so  no  other  empirical  relation  or  reference 
data  from  experiments  is  needed  in  the  numerical  calculations.  The  only 
parameter  which  needs  to  be  specified,  in  addition  to  the  heat  pipe 
dimensions  and  the  material  properties,  is  the  heat  transfer  coefficient  or 
the  radiation  emissivity  e  at  the  outer  wall  surface  of  the  condenser. 
This,  we  believe,  is  one  of  the  major  accomplishments  in  the  numerical 
analysis  of  heat  pipe  dynamics. 

6.6  C0NCLPSI0NS 

The  liquid- vapor  interface  energy  equation  (6.6)  is  derived  and  the 
conduction  model  for  the  heat  transfer  in  the  wick  proves  to  be  a 
reasonable  simplification  for  the  analysis  of  the  transient  behavior  of  a 
heat  pipe.  For  the  compressibility  of  the  vapor  flow,  the  pressure  is 
chosen  as  a  dependent  variable  and  the  equation  of  state  is  directly 
applied  during  the  iterations.  The  numerical  results  of  the  simulated 
vapor  flow  agree  well  with  the  corresponding  experimental  data.  The 
numerical  model  and  computer  code  can  properly  predict  the  transient 
process  of  a  heat  pipe.  The  transient  response  of  the  heat  pipe  to  a 
sudden  increase  of  the  heat  input  shows  that  the  wick  and  wall  play  a 
crucial  role  to  the  transient  process  of  the  heat  pipe  and  it  is  very 
important  to  study  the  heat  pipe  as  a  conjugate  problem. 
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Section  VII 


ANALYSIS  OF  THE  TRANSIENT  COMPRESSIBLE 
VAPOR  FLOW  IN  HEAT  PIPES 


7.1  SUMMARY 

The  transient  compressible  one- dimensional  vapor  flow  dynamics  in  a 
heat  pipe  is  modeled.  The  numerical  results  are  obtained  by  using  the 
implicit  non-iterative  Beam- Warming  finite  difference  method.  The  model  is 
tested  for  simulated  heat  pipe  vapor  flow  and  actual  vapor  flow  in 
cylindrical  heat  pipes.  A  good  comparison  of  the  present  transient  results 
for  the  simulated  heat  pipe  vapor  flow  with  the  previous  results  of  a 
two-dimensional  numerical  model  is  achieved  and  the  steady  state  results 
are  in  agreement  with  the  existing  experimental  data.  The  transient 
behavior  of  the  vapor  flow  under  subsonic,  sonic,  and  supersonic  speeds  and 
high  mass  flow  rates  are  successfully  predicted.  The  one- dimensional  model 
also  describes  the  vapor  flow  dynamics  in  cylindrical  heat  pipes  at  high 
temperatures. 
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7.2  imODUCTION 

In  a  heat  pipe,  the  change  of  phase  of  the  working  fluid  in  the  closed 
system  is  used  instead  of  a  large  temperature  gradient  to  transport  a  large 
amount  of  energy.  The  attention  of  many  scientists  has  focused  on  this 
unique  phenomenon  since  the  concept  was  introduced.  The  vapor  from  the 
evaporator  carries  energy  to  the  condenser,  so  the  vapor  flow  in  the  core 
region  of  the  heat  pipe  plays  an  important  role  in  transferring  energy  from 
source  to  sink.  ' 

Many  researchers  have  studied  the  steady  one-  dimensional  compressible  * 

(Levy,  1968;  Brovalsky  et  al.,  1976;  Faghri,  1989;  and  Jang,  1988)  and  the 
steady  two-dimensional  vapor  flow  in  heat  pipes  (Bankston  and  Smith,  1972; 

Tien  and  Rohani,  1974;  Ooijen  and  Hoogendoorn,  1979;  Faghri,  1986;  and 
Faghri  and  Parvani,  1988).  The  common  cross  sections  of  the  vapor  space 
are  circular,  rectangular  (Jang,  1988;  Ooijen  and  Hoogendoorn,  1979),  and 
annular  (Faghri,  1986;  Faghri,  1989;  and  Faghri  and  Parvani,  1988)  and  are 
chosen  based  on  the  particular  application.  The  heat  flux  distributions  on 
the  surface  of  the  evaporator  and  condenser  are  uniform  except  for  those 
presented  by  Jang  (1988).  Recently,  the  transient  two-dimensional 
compressible  simulated  vapor  flow  in  heat  pipes  (Bowman,  1987;  and  Bowman 
and  Hitchcock,  1988)  was  solved  numerically,  and  the  experimental  data  for  * 

the  steady  state  simulated  heat  pipe  vapor  flow  was  obtained  by  Bowman 
(1987).  However,  the  vapor  flow  was  not  actual  vapor  flow  in  a  heat  pipe  * 

and  the  numerical  and  experimental  data  were  presented  only  in  terms  of  the 
pressure.  The  transient  two-dimensional  compressible  vapor  flow  in  a  heat 
pipe  with  a  rectangular  cross  section  also  was  studied  numerically  by 
Issacci  et  al.  (1988),  but  only  the  axial  and  radial  velocity  profiles  were 
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presented.  No  comparison  with  the  existing  experimental  data  was  made  and 
the  transient  behavior  of  vapor  flow  in  a  heat  pipe  was  not  described. 

During  the  start-up  of  high  temperature  heat  pipes,  the  extremely 
small  density  of  the  vapor  causes  the  vapor  flow  to  attain  sonic  and 
supersonic  velocities  for  a  relatively  small  heat  input.  Thus,  the  correct 
description  of  the  transient  vapor  flow  is  essential  to  predict  the 
successful  start-up  and  to  estimate  the  overall  performance  of  the  entire 
heat  pipe.  The  governing  equations  for  the  vapor  flow  as  well  as  those  for 
the  wall  and  wick  regions  should  be  solved  simultaneously.  Also,  the 
development  of  the  one- dimensional  transient  model  for  the  vapor  flow  has 
been  suggested  due  to  the  large  amount  of  computer  time  required  for  the 
two-dimensional  model  (Bowman,  1987;  and  Bowman  and  Hitchcock,  1988).  For 
this  purpose,  a  simple  and  efficient  mathematical  model  is  desired  for  each 
region.  All  previous  one- dimensional  models  for  vapor  flow  are  for  the 
steady  state  condition  and  the  viscous  dissipation  was  neglected,  which  is 
important  for  high  temperature  applications. 

This  paper  describes  the  mathematical  model  and  the  numerical  method 
of  solution  for  the  transient  compressible  one- dimensional  vapor  flow  in 
the  heat  pipe.  A  comparison  of  the  numerical  results  with  the  simulated 
transient  two-dimensional  numerical  results  and  experimental  data  for  the 
steady  state  given  by  Bowman  (1987)  is  also  presented.  In  addition,  the 
numerical  results  from  the  present  model  for  the  actual  vapor  flow  in  the 
cylindrical  high  temperature  heat  pipe  are  compared  with  the  experimental 
data  obtained  at  the  steady  state  by  Ivanovskii  et  al.  (1982). 
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7.3  MATifEMATTCAL  MODELING 

The  one- dimensional  transient  compressible  vapor  flow  is  considered  to 

predict  the  vapor  flow  dynamics  in  the  heat  pipe.  Even  though  a  uniform 

velocity  is  used,  the  friction  at  the  interface  is  incorporated  by  using 

the  expressions  for  the  friction  coefficients  which  are  found  from  the 

two-dimensional  numerical  results  given  by  Bowman  (1987).  The  viscous 

dissipation  in  the  vapor  region  is  included  and  the  vapor  is  assumed  to  be 

a  perfect  gas.  The  governing  equations  for  the  vapor  flow  in  the  heat  pipe  1 

with  negligible  body  forces  are  formulated  by  using  the  principles  of  the 

conservation  of  mass,  momentum,  and  energy  in  a  control  volume  of 

o 

cross-sectional  area,  t  Dy  /4,  and  width,  dx. 

The  governing  equations  are  written  in  a  compact  vectorial  form  as 
follows: 


dd 

M  + 


(7.1) 


where 
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where  VQ(x)  refers  to  the  velocity  at  the  wall  with  a  positive  value  for 
injection  and  a  negative  value  for  suction. 


The  equation  of  state  is  employed  to  relate  the  density,  pressure,  and 
temperature  in  the  vapor  space  as  follows: 


p  =  pRT 


(7.6) 


For  the  simulated  heat  pipe  vapor  flow,  the  temperature  was  evaluated  by 
using  the  equation  of  state  because  a  change  of  phase  was  not  involved. 
For  the  actual  vapor  flow  in  the  cylindrical  heat  pipe,  the 
Clausius-  Clapeyron  relationship  was  used  to  predict  the  saturation 
temperature  of  the  vapor  from  the  pressure  as  given  by 


“in  E- 
pc 


(7.7) 


The  known  boundary  conditions  at  the  ends  of  the  heat  pipe  are  as 
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follows: 


U  =  0  at  x  =  0  and  L  (7.8) 

g  =  0  at  x  =  0  and  L  (7.9) 

The  conditions  for  the  density  and  pressure  at  the  ends  of  the  heat  pipe 
are  unknown,  so  physically  realistic  boundary  conditions  should  be  derived. 
In  general,  near  the  ends  of  the  heat  pipe  the  mass  flow  rate  is  small,  so 
the  axial  gradients  of  the  pressure  and  density  are  small.  In  the  region 
adjacent  to  the  exit  of  the  evaporator,  the  variations  of  the  pressure  and 
density  are  large.  Thus,  the  boundary  conditions  for  the  pressure  and 
density  can  be  assumed  as  follows: 


x  =  0  and  L 


(7.10) 


x  =  0  and  L 


(7.11) 


Bowman  (1987)  introduced  a  correlation  between  the  mass  flux  at  the 
wall  and  the  pressure  drop  across  a  porous  tube  wall  based  on  experimental 
measurements  for  the  simulated  heat  pipe  vapor  flow.  The  mass  flux  (p0^0) 
at  the  wall  for  the  simulated  heat  pipe  vapor  flow  was  evaluated  by  using 
this  correlation: 

A(p2)  =  3.639  x  109(p0Vo)2  ♦  1.7015  x  108(poVo)  (7.12) 

o 

where  A(p  )  is  the  absolute  value  of  the  difference  between  the  square  of 


the  uniform  source  pressure  and  the  square  of  the  vapor  pressure  in  the 

o 

blowing  section.  In  the  suction  section,  A(p  )  is  the  absolute  value  of 
the  difference  between  the  square  of  the  vapor  pressure  and  the  square  of 
the  uniform  sink  pressure.  A  change  of  phase  of  the  working  substance  was 
not  involved.  The  uniform  source  temperature  of  300  K  was  used  for  his 
experiment,  but  the  sink  temperature  was  not  specified.  To  evaluate  the 
terms  in  braces  {}  in  eqn.  (7.5),  the  source  temperature  is  used  at  the 
blowing  section  and  the  vapor  temperature  is  employed  at  the  suction 
section. 

Since  the  working  fluid  changes  phase  at  the  vapor- liquid  interface  in 
actual  heat  pipes,  the  temperature  at  the  interface  is  the  saturation 
temperature,  but  the  temperature  in  the  vapor  space  may  be  quite  different 
from  the  saturation  temperature  for  high  temperature  heat  pipes.  For  the 
one- dimensional  model,  the  properties  are  the  area- averaged  properties  so 
that  the  temperature  in  the  vapor  space  is  not  the  interface  temperature 
but  is  also  not  the  saturation  temperature.  The  vapor  temperature  can  be 
evaluated  from  the  energy  equation  and  the  saturation  temperature 
corresponding  the  vapor  pressure  can  be  obtained  from  the 

Clausius- Clapeyron  relationship.  However,  this  saturation  temperature  is 
not  the  actual  interface  temperature  either.  The  correct  estimation  of  the 
terms  in  braces  {}  in  eqn.  (7.5)  is  uncertain  due  to  using  the  area- 
*  averaged  properties.  Since  a  heat  pipe  is  a  closed  system,  the  application 

of  the  correct  values  of  heat  input  and  output  is  important.  To  eliminate 
this  difficulty,  terms  in  braces  {}  in  eqn.  (7.5)  are  replaced  by  using  the 
heat  flux  applied  on  the  surfaces  of  the  evaporator  and  condenser  as 
follows: 
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r  v  (x)  -] 

q  =  P0V0(x)  [  h0(x)  +  -J -  j  (7.13) 

7.4  FRICTION  COEFFICIENTS 

Bowman  (1987)  measured  the  turbulence  intensity  in  the  simulated  heat 
pipe  to  determine  the  characteristics  of  the  vapor  flow.  For  axial 
Reynolds  numbers  up  to  Re  =  10°,  laminar  flow  was  observed  in  the  blowing 
section  and  was  retained  in  the  suction  section  for  axial  Reynolds  numbers 
less  than  Re  -  12000.  The  transition  from  laminar  to  fully  turbulent  was 
predicted  in  the  entrance  region  of  the  suction  section  for  axial  Reynolds 
numbers  greater  than  Re  =  12000.  Unlike  flow  in  impermeable  tubes,  laminar 
flow  was  maintained  for  axial  Reynolds  numbers  greater  than  2000.  For  the 
supersonic  case  in  the  suction  and  blowing  section  the  flow  remained 
laminar  until  a  shock  wave  occurred  and  then  turbulent  flow  abruptly 
appeared,  which  showed  that  no  transition  region  existed. 

Since  the  mathematical  model  is  one- dimensional ,  proper  expressions 
for  the  friction  coefficient  are  necessary  to  take  into  account  the 
frictional  losses.  Bowman  (1987)  carried  out  eleven  numerical  simulations 
by  using  the  two-dimensional  numerical  model  to  evaluate  the  friction 
coefficients  according  to  the  characteristics  of  the  vapor  flow.  The 
following  expressions  for  the  friction  coefficients  including  the 
compressibility  effect  were  correlated  by  using  the  friction  coefficients 
obtained  from  the  two-dimensional  model.  The  friction  coefficient  for 
laminar  flow  in  the  condenser  or  evaporator  is 

f  =  (1.2337  -  0.2337e'0,0363Reo)e1'2M2  (7.14) 
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The  absolute  value  of  the  radial  Reynolds  number  at  the  wall,  ReQ,  is  used 
in  the  evaporator  and  condenser.  For  the  adiabatic  section  where  the  wall 
radial  Reynolds  number  is  zero,  eqn.  (7.14)  is  identical  to  that  for  the 
impermeable  circular  tube.  For  fully- developed  turbulent  flow  in  the 
condenser,  the  friction  coefficient  is 


,  _  0.046 
1  Re4'2 


1  +  55fte 


Vo« 


(7.15) 


For  transition  flow  in  the  condenser  entrance  region,  the  friction 
coefficient  is  defined  as 


-0.412x2 

e 


(7.16) 


For  the  one- dimensional  numerical  model  presented  in  this  paper,  the 
vapor  flow  in  the  evaporator  and  adiabatic  sections  is  assumed  to  be 
laminar.  In  the  condenser,  laminar  flow  is  assumed  for  the  axial  Reynolds 
numbers  below  Re  =  12000  at  the  entrance  of  the  condenser  and  for 

supersonic  flow.  When  the  axial  Reynolds  number  is  greater  than  Re  = 

12000,  transition  and  turbulent  flows  are  considered  in  the  condenser. 

Also,  after  a  shock  wave  turbulent  flow  is  assumed  in  the  condenser. 
Equation  (7.14)  is  used  to  evaluate  the  friction  coefficients  in  the 
\  evaporator  and  adiabatic  sections  and  is  also  employed  in  the  condenser 

when  the  axial  Reynolds  number  at  the  entrance  of  the  condenser  is  less 

than  Re  =  12000  and  the  vapor  flow  is  supersonic.  When  the  axial  Reynolds 
number  is  larger  than  Re  =  12000,  eqn.  (7.16)  is  applied  to  the  transition 
region,  which  is  assumed  to  exist  from  the  entrance  of  the  condenser  to 
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about  60  percent  of  the  condenser  length  based  on  the  experimental  data 
shown  in  Fig.  3.8  given  by  Bowman  (1987).  Equation  (7.15)  is  used  for  the 
turbulent  flow  in  the  condenser.  The  initial  location  of  the  transition 
region,  x.  . ,  is  equal  to  the  location  at  the  entrance  of  the  condenser, 
and  from  Fig.  7.3.8  given  by  Bowman  (1987),  Xp  _  and  Xp  _  ^  were 
estimated  to  be  0.7  and  0.6,  respectively. 

7.5  NUMERICAL  FORMULATION 

From  the  many  schemes  (Anderson  et  al.,  1984)  available  for  the 
solution  of  the  compressible  flow  problem,  the  Beam- Warming  finite 
difference  scheme  is  chosen  to  transform  the  governing  equation  (7.1)  to 
the  finite  difference  formulation.  This  scheme  is  a  non-iterative  implicit 
method  and  is  similar  to  ADI  for  multidimensional  flow  problems  by  using 
the  factorization  which  retains  the  tridiagonal  block  matrix. 


The  spatial  derivatives  are  approximated  by  using  the  three- point 
second- order  accurate  central  difference  approximation  for  the  interior 
points  and  the  one-sided  second- order  accurate  difference  approximation  for 
the  boundary  nodes.  After  the  approximation  operators  are  applied,  the 
system  of  equations  becomes 


[RHS]ni  i 


2  I 
z  ’  ’  max- 1 

(7.17) 


where 

tfnD  =  Dn+1  -  Dn 
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where  [I]  is  the  unit  matrix,  [A] ,  [P] ,  [R]  and  [S]  are  the  Jacobian 

matrices,  and  c  is  the  coefficient  of  the  dissipative  term  to  damp  the 
oscillation. 


This  time  difference  formula  reproduces  many  different  schemes  with 
the  appropriate  choice  of  9^  and  The  scheme  is  second- order  accurate  in 
time  when  9^  =  1/2  +  9 2  and  first- order  accurate  otherwise.  For  0^-1  and 
#2  =  1/2,  the  formula  becomes  second- order  accurate  in  time  over  three  grid 
points.  The  system  of  equations,  (7.17),  has  the  following  block 
tridiagonal  structure: 

[JKL]  {<JnD.}  =  {RHSin}  (7.18) 

where  [JKL]  represents  the  banded  coefficient  matrix  of  which  components 
are  3x3  matrices  for  one- dimensional  vapor  flow,  and  {<5nD^}  and  {RHS^n} 
are  column  vectors.  The  tridiagonal  block  matrix  size  is  now  (3  x  Imax  2) 
x  (3  x  I  0)  where  I  is  the  number  of  nodal  points.  This  system  of 
equations  can  be  solved  using  the  conventional  methods  for  solving  block 
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tridiagonal  systems  of  equations.  The  vector  of  unknowns  at  n  +  1  time 
step  is  then  determined  by  simply  adding  £nD  to  the  value  of  Dn  at  n  step. 
The  primitive  variables  ( p ,  U,  p,  T)  can  be  obtained  from  Dn+*. 

A  total  of  80  nodes,  which  is  the  minimum  number  of  nodes  to  obtain 
accurate  results,  are  evenly  spaced  in  the  axial  direction  and  a  time  step 
of  0.1  x  10  s  is  used  for  the  simulated  heat  pipe  vapor  flow.  Since  the 
heat  flux  at  the  wall  in  the  evaporator  is  different  from  that  in  the 
condenser  due  to  the  different  lengths  for  the  cylindrical  heat  pipe,  the 
coarse  nodal  system  presented  some  difficulty  to  reach  the  steady  state. 
For  the  cylindrical  heat  pipe,  200  evenly- spaced  nodes  are  used  in  the 

_  Q 

axial  direction  and  a  time  step  of  0.1  x  10  s  is  employed. 

7.6  RESULTS  AWT)  DISCUSSION 

7.6.1  Comparison  with  the  simulated  heat  pipe  vapor  flow 
A  comparison  of  the  numerical  results  with  the  experimental  data  given 
by  Bowman  (1987)  is  desired  to  verify  the  mathematical  model  and  algorithm. 
However,  the  existing  experimental  data  was  obtained  by  simulating  the 
vapor  flow  of  a  cylindrical  heat  pipe  with  a  porous  pipe  which  has  an 
inside  diameter  of  1.65  cm  and  a  length  of  0.61  m  as  shown  in  Fig.  7.1. 
The  blowing  and  suction  sections  have  equal  lengths  and  were  simulated  by 
the  injection  and  suction  of  air  without  phase  change  at  the  interface. 
Also,  uniform  source  and  sink  pressures  for  the  blowing  and  suction  regions 
are  specified  instead  of  the  radial  mass  flow  rate.  Thus,  to  simulate  the 
experiment  eqn.  (7.12)  is  used  to  obtain  the  radial  mass  flux  with  the 
known  source  and  sink  pressures. 
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AIR  IN 


Figure  7.1.  Schematic  diagram  of  model  for  air  flow  in  the  porous  pipe 
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7.6. 1.1  Transient  results 

The  existing  experimental  data  was  obtained  at  the  steady  state  so  the 

present  transient  numerical  results  can  only  be  compared  with  the  transient 

numerical  results  for  the  two-dimensional  model  (Bowman,  1987).  For  this 

purpose,  the  same  geometry  and  physical  conditions  are  used  such  as  the 

5  2 

source  pressure  of  2.06  x  10  N/m  (30  psia)  and  the  sink  pressure  of  1.03 
5  2 

x  10  N/m  (15  psia)  corresponding  to  case  B.l. 

Initially,  the  velocity  of  the  vapor  is  zero  and  the  pressure  and 
temperature  are  the  same  as  the  source  pressure  and  temperature, 
respectively.  To  simulate  the  transient  flow,  the  sink  pressure  is 
suddenly  lowered  to  1.03  x  10  N/m  (15  psia),  while  the  source  pressure 
remains  the  same  as  the  initial  pressure.  This  difference  between  the 
source  and  sink  pressures  initiates  the  air  flow  from  the  blowing  section 
to  the  suction  section.  The  pressures  at  the  inlet  of  the  blowing  section, 
the  center  of  the  pipe,  and  the  end  of  the  suction  section  are  plotted  to 
compare  with  the  numerical  results  for  the  two-dimensional  model  as  shown 
in  Fig.  7.2. 

Figure  7.2  shows  the  following  transient  behavior  of  the  air  flow  in 
the  porous  pipe.  Since  the  sink  pressure  is  abruptly  changed  from  the 
initial  pressure  to  1.03  x  10  N/m  (15  psia)  along  the  entire  suction 
section,  the  pressures  at  the  center  and  last  nodes  decrease  immediately 
due  to  the  evacuation  of  air.  For  this  period,  the  blowing  section 
pressures  adjacent  to  the  suction  section  start  to  decrease  due  to  the  flow 
of  mass  from  the  blowing  section  to  the  suction  section,  but  the  pressure 
near  the  beginning  of  the  blowing  section  remains  constant.  Also,  the  mass 
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Present  1-D  transient  model 
for  Case  B.l 


Figure  7.2.  Comparison  of  the  present  numerical  results  with 
Bowman’s  2-d  model  for  pressure  variations  with  time  at 
three  locations  of  the  porous  pipe:  Case  B.l 


flow  rate  from  the  blowing  section  to  the  suction  section  is  not  sufficient 
to  influence  the  end  of  the  suction  section  so  that  the  pressure  at  this 
point  decreases  faster  than  that  at  the  center  of  pipe. 

At  about  0.8  -  1.0  x  10  s,  the  pressure  at  the  end  of  the  suction 
section  reaches  the  minimum  value  and  then  starts  to  increase  while  the 
pressure  at  the  center  of  pipe  keeps  decreasing  due  to  the  frictional  less 
and  the  acceleration  of  the  flow.  At  this  time,  the  pressures  over  the 
entire  blowing  section  become  less  than  the  initial  pressure  so  that  the 
mass  flow  rate  is  sufficient  to  influence  the  end  of  suction  section.  As 
the  pressure  in  the  blowing  section  decreases  and  the  source  pressure 
remains  constant,  the  mass  flow  rate  from  the  blowing  section  to  the 
suction  section  increases.  Thus,  the  pressure  at  the  center  node  keeps 
decreasing  and  the  pressure  at  the  end  of  the  suction  section  rises  due  to 
the  contribution  of  the  mass  from  the  blowing  section. 

.  3 

At  about  3.5  x  10  s,  the  pressures  at  all  three  points  reach  the 
steady  state.  As  expected,  the  pressure  at  the  end  of  the  suction  section 
does  not  recover  completely  due  to  the  frictional  loss  at  the  pipe  wall. 
Figure  7.2  shows  that  the  present  results  and  numerical  results  for  the 
two-dimensional  model  (Bowman,  1987)  are  in  agreement. 

7. 6. 1.2  Steady  state  results 

When  the  present  numerical  results  reach  the  steady  state,  those 
results  are  compared  with  the  experimental  data  given  by  Bowman  (1987) . 
Numerical  calculations  are  conducted  for  four  different  sets  (i.e.,  cases 
B.2,  B.3,  B.4,  and  B.5)  of  the  source  and  sink  pressures  by  using  eqns. 


(7.14  -  7.16)  for  the  friction  coefficient.  The  cases  examined  are  as 
follows: 

Case  B.2:  PgQ  =  1.39  x  105  and  Pgk  =  1.21  x  105  N/m2 
Case  B.3:  P  -  2.09  x  105  and  Pgk  =  1.57  x  105  N/m2 
Case  B.4:  P  =  2.68  x  105  and  Pgk  =  1.69  x  105  N/m2 
Case  B.5:  Pg0  =  5.30  x  105  and  Pgk  =  1.06  x  105  N/m2 

Figure  7.3  shows  the  comparison  of  the  pressure  distributions  along  the 
axial  direction.  The  top  three  lines  in  Fig.  7.3  show  the  pressure 
distributions  for  the  low  mass  flow  rates.  The  pressures  in  the  blowing 
section  decrease  due  to  friction  and  the  acceleration  of  the  flow  caused  by 
mass  injection,  but  the  pressures  in  the  suction  section  increase  owing  to 
the  deceleration  of  the  flow  by  the  extraction  of  mass.  However,  the 
pressures  at  the  end  of  the  suction  section  are  less  than  those  at  the 
beginning  of  the  blowing  section  because  of  the  loss  due  to  friction.  In 
these  three  cases,  the  pressure  distributions  at  the  steady  state 
correspond  well  to  those  for  low  temperature  heat  pipes  and  are  in 
agreement  with  the  experimental  data. 

The  fourth  data  profile  (case  B.5)  shows  the  pressure  variation  in  the 
axial  direction  for  the  high  mass  flow  rate.  Unlike  the  previous  three 
cases,  the  pressure  drop  in  the  blowing  section  is  very  large.  The 
pressure  ratio  at  the  exit  of  the  blowing  section  is  about  0.4  and  this 
ratio  corresponds  to  a  Mach  number  of  M  =  1.  After  the  pressure  decreases 
in  the  blowing  section,  the  pressure  keeps  decreasing  in  the  entrance 
region  of  the  suction  section  due  to  the  expansion  of  air  even  though  mass 
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o  =  Case  B.2 
•  =  Case  B.3 


7.3.  Comparison  of  the  present  numerical  results  with 
experimental  pressure  variations  (Bowman,  1987) 
in  the  porous  pipe 


removal  occurs.  Then,  the  pressure  suddenly  increases  and  then  continues 
to  increase  as  the  flow  slows  down.  This  implies  that  a  shock  wave  occurs 
at  the  place  where  the  pressure  changes  abruptly.  When  a  shock  wave  does 
not  exist  in  the  suction  section,  the  pressure  is  supposed  to  decrease 
along  the  suction  section.  The  one- dimensional  model  predicts  the 
supersonic  flow  and  shock  wave  in  the  suction  section  and  the  comparison  of 
the  numerical  results  and  experimental  data  shows  a  good  agreement  except 
^  for  the  region  immediately  after  the  shock  wave. 

V  The  variations  of  the  pressure,  temperature,  velocity  and  density  at 

the  steady  state  corresponding  to  cases  B.4  and  B.5  are  shown  in  Figs.  7.4 
and  7.5,  respectively.  For  the  low  mass  flow  rate,  the  temperature  and 
density  in  the  blowing  section  decrease  corresponding  to  the  decrease  in 
pressure  and  the  velocity  increases  due  to  the  mass  injection  as  shown  in 
Fig.  7.4.  The  Mach  number,  however,  is  less  than  M  =  1  at  the  exit  of  the 
blowing  section,  so  the  velocity  in  the  suction  section  decreases  because 
of  the  extraction  of  mass.  Also,  the  temperature  and  density  increase  in 
the  suction  section.  Figure  7.5  shows  the  axial  variations  of  the 
temperature,  pressure,  velocity  and  density  for  the  high  mass  flow  rate. 
After  the  sonic  velocity  is  reached  at  the  exit  of  the  blowing  section,  the 
)  velocity  keeps  increasing  in  the  suction  section  until  a  shock  wave  occurs. 

Then,  the  velocity  decreases  to  U  =  0  at  the  end  of  the  pipe.  As  shown  by 
}  the  density  profile,  the  air  expands  near  the  entrance  region  of  the 

suction  section  and  then  the  density  suddenly  increases  after  the  shock 
wave. 


The  most  interesting  aspect  in  Figs.  7.4  and  7.5  is  the  variation  of 
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Figure  7.4.  Axial  variations  of  temperature,  pressure,  velocity, 

and  density  for  Case  B.4 
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Figure  7.5.  Axial  variations  of  temperature,  pressure,  velocity, 

and  density  for  Case  B.5 


temperature.  The  temperature  at  the  end  of  the  suction  section  is  greater 
than  that  at  the  beginning  of  the  blowing  section.  In  the  present  model 
viscous  dissipation  is  included.  The  Mach  number  at  the  entrance  of 
suction  section  is  about  M=0.5  for  case  B.4.  Since  the  present  model  is 
one- dimensional ,  the  derivative  of  the  axial  velocity  with  respect  to  the 
radius  is  zero,  but  the  friction  effect  at  the  interface  between  the  wall 
and  the  vapor  is  included  by  using  eqns.  (7.14  -  7.16)  for  the  friction 
factor.  This  effect  corresponds  to  viscous  dissipation  due  to  the  axial 
velocity  derivative  with  respect  to  the  radius.  The  increase  in 
temperature  is  due  to  viscous  dissipation. 

7.6.2  Comparison  with  actual  vapor  flow  in  cylindrical  heat  pipes 

The  model  was  tested  for  the  actual  vapor  flow  in  a  sodium  heat  pipe 
corresponding  to  the  experiment  given  by  Ivanovskii  et  al.  (1982)  as  shown 
in  Fig.  7.6.  The  lengths  ofthe  evaporator,  adiabatic,  and  condenser 
sections  are  0.1,  0.05,  and  0.35  m,  respectively.  The  diameter  of  the 
vapor  space  is  0.014  m.  The  experimental  data  (Ivanovskii  et  al.,  1982) 
represents  only  the  saturation  temperature  and  heat  transfer  rate  (Q  =  560 
W)  at  the  steady  state.  The  initial  conditions  and  actual  boundary 
conditions  applied  on  the  surface  of  the  evaporator  and  condenser  are 
unknown.  The  present  numerical  model  is  transient,  however,  so  the  initial 
conditions  and  history  of  heat  input  and  output  in  the  evaporator  and 
condenser  sections  are  needed. 

For  the  numerical  calculations,  the  uniform  initial  temperature  of  810 
K  is  used  and  the  vapor  is  assumed  to  be  saturated  at  the  initial 
temperature.  Initially,  the  velocity  of  the  vapor  is  zero.  The  same 
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Lc  =  0.35  ra 


Figure  7.6.  Schematic  diagram  of  model  for  sodium  vapor  flow 

in  the  heat  pipe 


amount  of  heat  (Q  =  560  ¥)  is  uniformly  applied  on  the  surface  of  the 
evaporator  and  the  convective  boundary  condition  is  '^ed  on  the  surface  of 
the  condenser.  The  reference  temperature  of  300  K  for  the  convective 
boundary  condition  is  employed  and  the  heat  transfer  coefficient  is 
determined  by  iteration.  At  first,  an  arbitrary  initial  heat  transfer 
coefficient  is  guessed.  Vhen  the  numerical  results  reach  the  steady  state, 
the  saturation  temperature  at  the  end  cap  of  the  evaporator  is  compared 
with  the  experimental  data  at  the  same  location.  This  procedure  is 
repeated  until  the  same  temperature  is  obtained  at  the  steady  state.  For 
this  test,  the  heat  transfer  coefficient  of  69.1  W/m  -  K  is  used. 

Figure  7.7  shows  the  axial  variation  of  the  saturation  and  vapor 
temperatures,  pressure,  velocity  and  density  obtained  from  the  present 
numerical  model  and  the  experimentally  measured  saturation  temperature 
distribution  (Ivanovskii  et  al.,  1982).  The  pressure,  temperature  and 
density  in  the  evaporator  decrease  and  the  velocity  increases  due  to  the 
injection  of  mass  and  the  effect  of  friction.  In  the  adiabatic  section, 
the  vapor  temperature  increases  because  of  viscous  dissipation  and  the 
pressure  decreases  owing  to  friction  at  the  interface.  The  density  also 
decreases  while  the  velocity  continues  to  increase.  The  vapor  temperature 
in  the  condenser  continues  to  increase  due  to  viscous  dissipation  so  that 
the  vapor  temperature  at  the  end  cap  of  the  evaporator  is  less  than  that  at 
the  end  cap  of  the  condenser.  The  pressure  recovery  in  the  condenser  is 
almost  negligible.  This  may  result  from  dominant  friction  effect  at  the 
interface  of  the  condenser  compared  to  the  effect  of  mass  extraction  in 
this  long  condenser.  In  the  adiabatic  section  a  difference  between  the 
calculated  and  measured  saturation  temperatures  is  observed,  but  the  trend 
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Figure  7.7.  Axial  variations  of  temperature,  pressure,  velocity, 
and  density  of  the  sodium  heat  pipe  at  steady  state 


of  the  saturation  temperature  variation  is  the  same.  The  trend  of  the 
vapor  temperature  in  the  condenser  is  quite  different  fiuiu  that  of  the 
saturation  temperature  in  the  same  region,  so  the  saturation  temperature 
may  not  be  assumed  to  be  the  vapor  temperature  in  the  condenser  for  the 
one- dimensional  model. 

The  effect  of  viscous  dissipation  in  the  vapor  flow  for  the  high 
temperature  heat  pipe  is  investigated  as  shown  in  Fig.  7.8.  The  vapor  ' 

temperature  at  the  end  cap  of  the  evaporator  with  viscous  dissipation  is 
less  than  that  at  the  end  cap  of  the  condenser.  However,  the  vapor  ) 

temperature  at  the  end  cap  of  the  evaporator  without  viscous  dissipation  is 
almost  the  same  as  that  at  the  end  cap  of  the  condenser.  In  the  adiabatic 
section,  the  vapor  temperature  with  viscous  dissipation  increases  while 
that  without  viscous  dissipation  does  not  change.  Therefore,  the  viscous 
dissipation  should  be  taken  into  account  in  the  energy  equation  for  the 
high  temperature  heat  pipe. 

7.7  CONCLUSIONS 

A  model  for  the  transient  one- dimensional  compressible  vapor  flow  in 
the  cylindrical  heat  pipe  is  developed.  This  model  predicts  the  vapor  flow 
in  cylindrical  heat  pipes  as  well  as  simulated  heat  pipes  for  the  subsonic,  V 

sonic,  and  supersonic  flows  under  transient  and  steady  state  conditions. 

The  vapor  flow  quickly  reaches  the  steady  state  condition.  The  4 

distributions  of  the  temperature  and  pressure  during  the  transient  state 
are  quite  different  from  those  for  the  steady  state.  The  viscous 
dissipation  terms  play  an  important  role  in  the  energy  equation  and  have  to 
be  taken  into  account.  The  one- dimensional  compressible  model  predicts  the 
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Figure  7.8.  Comparison  of  the  temperature  variations  with  an 
without  viscous  dissipation  for  actual  vapor  flow 


experimental  data  well  for  the  cylindrical  heat  pipe  and  the  simulated  heat 
pipe  at  the  steady  state.  The  experimental  data  for  the  transient  state 
are  needed  to  understand  clearly  the  transient  behavior  of  the  vapor  flow 
in  the  heat  pipe  both  at  low  and  high  temperatures.  The  one- dimensional 
model  can  reduce  the  computational  effort  needed  to  solve  the  vapor  flow 
problem. 
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